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ABSTRACT 


A FORTRAN IV computer program which incorporates the method of characteristics 
was written to assist in the design of supersonic inlets. The computer program was 
written with the intention of studying many types of inlet configurations with a minimum 
of computer time. Particular attention was given to a discussion of a reformulation of 
the boundary value problem to introduce throat Mach number and flow angle as direct 
program input quantities. Comparison between the computer results and experimental 
data indicates good agreement for a number of different configurations. 



DESIGN OF SUPERSONIC INLETS BY A COMPUTER PROGRAM 
INCORPORATING THE METHOD OF CHARACTERISTICS 
by Bernhard H. Anderson 
Lewis Research Center 

SUMMARY 

A FORTRAN IV computer program which incorporates the method of characteristics 
was written to assist in the design of supersonic inlets. The program is versatile so 
that many inlet types can be studied with a minimum of computer time. 

Particular attention was given to reformulating the boundary value problem so that 
the engineering design parameters, that is, throat Mach number and flow angle, could be 
introduced as direct input quantities to the computer program. Attention was also given 
to a discussion of streamline tracing by a numerical integration of the mass flux from a 
known physical boundary to determine the opposite boundary. 

Selected examples are presented and compared with experimental results to illus- 
trate and corroborate the numerical methods used. The report also includes FORTRAN 
listings for the main programs and subroutines and discussions of their various func- 
tions. 


INTRODUCTION 

The mathematical methods for designing supersonic inlets by characteristic theory 
have been available for many years (ref. 1). However, the major obstacle in the design 
of supersonic inlets by the method of characteristics has been the great complexity of 
the numerical procedures required for configurations other than simple wave solutions. 
Thus, to facilitate the design of inlets, computational programs have been written for 
use on high-speed computers (refs. 2 to 4). Because of the often complex nature of the 
boundaries which constitute a supersonic inlet, these programs have been limited to 
specific configurations. In addition, computational methods currently being used for the 
design of high-performance inlets often necessitate extensive trial-and-error procedures 
to arrive at an optimum design. This results from the fact that some of the primary 


input data to the computer program are mathematical in significance and only indirectly 
related to important engineering design parameters, such as throat Mach number, flow 
angle, and distortion. 

For these reasons, a FORTRAN IV computer program, employing the method of 
characteristics, was written with two objectives in mind: (1) to study a greater variety 
of supersonic inlet configurations and (2) to reduce the time required for trial -and-error 
procedures to arrive at optimum inlet design. The computer program was written with 
the intention of being able to construct a variety of inlet configurations by interchanging 
specific subroutines. In this manner, greater flexibility of choice was attained, and the 
time required to program a specific inlet configuration was greatly reduced. The second 
objective was accomplished by a reformulation of the boundary value problem for hyper- 
bolic equations. By this reformulation of the boundary data, the engineering design 
quantities, such as throat Mach number and flow angle, were introduced as direct input 
quantities to the computer program. As a consequence of introducing the engineering 
parameters as input, the computer program will calculate the surface contours required 
to satisfy the specified throat conditions with the minimum inviscid throat distortion. 

To facilitate the use of the method of characteristics in the design of supersonic in- 
lets, numerical procedures have been derived from the general theory and adapted to 
high-speed computers. Particular attention was given in this study to a discussion of 
boundary data and ways they may be applied for more effective inlet design. In addition, 
attention was given to a discussion of streamline tracing by a numerical integration of 
the mass flux from a known physical boundary. This technique was used extensively in 
the computer program with great effectiveness. 

To illustrate and corroborate the numerical methods, selected examples are pre- 
sented and compared with experimental results. Sample problems are also included to 
illustrate the use of the program to achieve specific objectives. This report also pre- 
sents FORTRAN listings for the main programs and subroutines, descriptions and func- 
tions of the subroutines, explanations of the input-output quantities, and general program 
procedures. 


SYMBOLS 

A dimensionless coefficient defined by eq. (3) 

B dimensionless coefficient defined by eq. (4) 

C characteristic 

C, additive drag coefficient 

CXy ci 

c dimensionless speed of sound (see appendix A) referenced to critical speed 
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D dimensionless coefficient defined by eq. (5) 

E dimensionless coefficient defined by eq. (6) 

F dimensionless coefficient defined by eq. (7) 

G dimensionless coefficient defined by eq. (7) 

M Mach number 

m mass flow 

P total pressure 

p static pressure 

q dimensionless local velocity relative to critical speed 

R coefficients for equation of boundary 

r dimensionless radius, + y^ 

u dimensionless x-component of velocity (see appendix A) relative to critical 

speed 

v dimensionless y-component of velocity (see appendix A) relative to critical 

speed 

x dimensionless x-coordinate relative to cowl lip radius 

y dimensionless y-coordinate relative to cowl lip radius 

a(x, y) parameter defining C + characteristic family (see eq. (A10)) 

|3(x,y) parameter defining C_ characteristic family (see eq. (A10)) 

0_ shock angle 

y ratio of specific heats 

6 deflection angle 

e tangent of solid boundary 

77 cotangent of conic ray half angle (see eq. (Bl)) 

X Eigen values of eq. (A7) 

0 conic ray half angle 

M local Mach angle with respect to local flow angle 

p dimensionless density relative to critical conditions 

0 local flow angle with respect to coordinate system 

? characteristic directions (see eq. (A15)) 
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\p stream function 

<j2 . parameter defined as (cr - 1)/ (cr + 1) 

r parameter defined by eq. (A6) 

Subscripts: 

c cowl surface 

cb centerbody surface 

ref reference conditions 

w wall conditions 

0 free -stream conditions 

Superscript: 

— average conditions 


GOVERNING EQUATIONS 

To facilitate the use of the characteristic theory in the design of supersonic inlets, 
numerical procedures were developed with the following underlying assumptions: 

(1) Steady three-dimensional flow with cylindrical symmetry 

(2) Homentropic flow within shockless flow regions 

Assumption (1) is not overly restrictive since a large class of two- and three- 
dimensional inlets of practical interest are satisfied by this condition; assumption (2) 
is implemented along physical boundaries by basing the calculations on the total- 
pressure recovery along that streamline. In general, the accuracy of results with 
assumption (2) was not appreciably different from isentropic calculations. In general, 
these errors were of the order of 0. 10 to 0. 50 percent in local Mach number. 


Characteristic Equations 

Flow with cylindrical symmetry is characterized by the conditions that all pertinent 
quantities depend only on the abscissa x along the axis and the distance y from the 
axis. The characteristic equations and compatibility conditions are found in reference 1 
and derived in appendix A. They are given, respectively, by 

dy = Adx 

dy = Bdx 
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where 


du = Ddv + Fdx C + > | 


du = Edv + Gdx C _ J 

. (2) 

A = tan (d + ju) 

(3) 

B = tan (S - IJ-) 

(4) 

D = -tan (5 - M) 

(5) 

E = -tan (3 + m) 

(6) 

F = G= ° 2 v 

(7) 


2 2 y 

u - c 3 


In addition to the above equations, Bernoulli’s Law gives the following relation: 


oV + v 2 ) + (1 - cr 2 )c 2 = 1 


( 8 ) 


where 


y + 1 


(9) 


The set of equations (1), (2), and (8) have been made dimensionless by dividing the spatial 
coordinates by the cowl lip radius and the velocity coordinates by the critical speed. The 
equations for two-dimensional flow are obtained by setting F = G = 0. 


Conical Flow Field Equations 

Further mathematical simplifications are possible for the case of axisymmetric flow 
if the flow field considered is that produced by a cone with sufficiently small half angle 
that the shock front is attached. Under these conditions, the velocity components u, v 
and speed of sound c depend only on the angle between the axis of symmetry and the 
family of concentric cones between the obstacle cone and shock wave. When the velocity 
components u and v are introduced as the independent and dependent variables (see 


5 


appendix B), respectively, the governing equation becomes 


vv 


uu 



(1 - or 2 ) (u + vv u ) 2 

i 2/2 2x 

1 - a (u + v ) 


( 10 ) 


where the subscripts indicate differentiation. Along the cone surface, the boundary con- 
dition is specified by the relation 


u 



(ID 


The conditions to be satisfied by the velocity immediately behind the conical shock are 
specified by the oblique shock equations 

u = (1 - a 2 )q n cos 2 P + — 
v = (qQ - u) cot 0 S 


u 

where q^ is the dimensionless free-stream velocity and /3 g is the shock -wave angle. 
Equations (10) to (14) are derived in reference 1 and listed in appendix B. 


( 12 ) 

(13) 

(14) 


Oblique Shock-Wave Relations 

The standard oblique shock relation used was obtained from reference 5 and is given 
by 


where 


sin® /3 + b sin^ P + c sin 2 j3 + d = 0 


b = - M — t. 2 - y sin 2 6 
M 2 


(15) 


( 16 ) 


6 


.( 17 ) 


c = 2“1li + 

M 4 


(y + 1\ 2 + y - 1 

V 2 / „2 


M J 


2 

sin 6 


d = - 


o 

cos 6 
M 4 


(18) 


M is the Mach number upstream of the shock wave, and 6 is the deflection angle. 


FORMULATION OF BOUNDARY DATA 

The type of data which must be prescribed to have a well-set problem is fundamental 
to establishing a numerical solution to the characteristic equations. In the design of 
supersonic inlets, a common situation which arises includes a known physical boundary 
along with an initial datum line which can be either noncharacteristic (such as a shock 
wave) or one member of the two-family set of characteristics (fig. 1). The solution of 
the flow field can then be established in the triangular domain ABC in the x,y plane. 

A more general way of formulating boundary requirements is to consider the flow 
field to be constrained between two boundaries (fig. 2), one of which may not be defined. 

In order to establish a well-set problem, data must be prescribed along the initial datum 
line and a pair of boundary data must be provided. The possible combinations of boundary 
data include 

(1) Prescribing the two boundary curves, y = y c ^(x) and y = y c (x), and solving for 

the flow field in the domain ABCD and conditions along each of the bounding 
surfaces 

(2) Prescribing one boundary surface, y = y c b(x), and the corresponding velocity 

distribution along the same surface, q = q c ^(x) (Under these conditions, the 
boundary value problem is to establish the flow field and the necessary contour 
of the second surface, y = y c (x), which are consistent with the initial data and 
prescribed boundary information. ) 

(3) Prescribing the velocity distribution for each of the surfaces, for example, 

q = q c j J (x) and q = q c (x) (Under these conditions, only the magnitude of the 
velocity is known on each of the two boundaries. The solution must thus include 
establishing the unknown surfaces and flow field which satisfy the prescribed 
initial and boundary data. ) 

For the purposes of this report, the three types of boundary data described are 
termed 

(1) Fixed -boundary problem (Case 1) 
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(2) Free-boundary problem (Case 2) 

(3) Doubly -free -boundary problem (Case 3) 

Fixed-Boundary Problem 

When both boundaries are fixed, the equations of the bounding surfaces become the 
controlling influence (i. e. , they are program input variables). For this situation, if a 
particular flow field must be established within the region ABCD (i. e. , throat Mach 
number and flow angle (fig. 2)), a trial -and-error iteration for the surfaces y = y cb (x) 
and y = y_(x) must be used. 

V 


Free-Boundary Problem 

The free-boundary problem is a useful way of prescribing boundary data if a par- 
ticular flow field is to be established, for example, in the throat region of a supersonic 
inlet. Prescribing the centerbody contour and velocity distribution establishes the 
manner in which the flow is to be compressed and the coordinates of the upper surface, 
y = y c (x), which is the streamline which passes through the cowl lip. To establish a 
uniform flow field downstream of domain ABCD, the velocity and surface angle of bound- 
ary (1) are held constant downstream from point C. This establishes the requirement of 
a particular throat Mach number with minimum inviscid flow distortions. 

For the free-boundary problem, the establishment of the unknown surface, y = y c (x), 
essentially reduces to finding the streamline which passes through point A in figure 2. 

In any inviscid flow, the streamlines, by definition, form surfaces across which there is 
no flow and consequently may be replaced by a solid boundary. The streamlines can be 
determined either by constructing them incrementally from known local conditions or by 
integration of the mass flux from known boundary conditions. This latter method was 
used extensively throughout the computer programs. Equations used for the streamline 
tracing are derived in appendixes C and D, and discussions as to their use appear in 
detail in later sections. 

The boundary for which the pair of data, y = y w (x) and q = q w (x), can be pre- 
scribed must satisfy the condition that the downstream leading characteristics from the 
initial datum line do not intersect that surface. If the wrong boundary is chosen, an 
overspecified problem will result. 
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Doubly-Free-Boundary Problem 


The doubly -free -boundary problem is an alternate way of prescribing boundary data 
and may have useful applications. Its disadvantage lies in the loss of direct control of 
flow direction, which becomes important in inlet design. This immediate disadvantage 
does not occur for either the fixed or free-boundary problems. In any case, the doubly- 
free -boundary problem was not explored further in this report. 


NUMERICAL PROCEDURES 

Each net point in the x,y plane is determined by the characteristic equations and 
compatibility conditions specified by equations (1) and (2). The properties u,v are 
obtained from the compatibility conditions, while the characteristic equations locate the 
point in the x,y plane. The computational work divides itself into four distinct groups; 
namely, 

(1) Conical flow field calculations 

(2) Basic characteristic field point 

(3) Basic characteristic boundary point 

(4) Shock point calculations 


Conical Flow Field Calculations 

The mathematical construction for the flow past an infinite cone with attached shock 

fronts is obtained from a numerical solution of equation (10) by means of a Runge-Kutta 

integration. The conical shock angle /3 is established by an iterative process which 

s 

begins with an initial estimate of the shock angle for the given free-stream conditions. 
The relations governing the transitions through conical shocks are the same as those for 
plane oblique shocks; the curvature of the shock cone has no effect. The necessary 
boundary conditions to be satisfied along the shock front can thus be described by equa- 
tions (12) to (14). A solution of equation (10) is found which is compatible with the shock 
values of u, v defined by equations (12) and (13) and has an initial slope specified by 
equation (14). The solution is then continued so that v u increases until 


along the cone surface. The cone angle thus obtained is compared to the specified cone 
angle, and new estimates of the shock angle 0 S are made until the cone angle and flow 
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angle are equal. When convergence has been reached, the boundary conditions u,v on 
the cone have been established. The solution for any point in the conical flow field can 
thus be calculated from equation (10) by specifying the ray angle associated with the flow 
field point. The formulae used for additive drag and mass flow spillage calculation in 
the conical flow field appear in appendix E. 


Basic Characteristic Field Point 

The program specifies each field point by a (K, I) index combination. Increasing 
the K -index (with I held constant) occurs along a C + characteristic (fig. 3), while in- 
creasing the I -index (with K held constant) indicates a C_ characteristic. The unknown 
position coordinates x(K, I), y(K, I) are obtained by a simultaneous solution of char- 
acteristic equations (1) in finite difference form, that is, 

x(K, I) = y(K, I - 1) - y(K - 1, I)j- Ax(K - 1, I) - Bx(K, I - 1) (19) 

A - B 


y(K, I) = y(K - 1, I) + A [x(K, I) - x(K - 1, I)] 


where 


A = 0. 50 [A(K, I) + A(K - 1, I)] 

B = 0. 50 [B(K, I) + B(K, I - 1)] 

The unknown velocity coordinates are obtained from a simultaneous solution of the com- 
patibility conditions (eq. (2)) in finite difference form, that is. 


v(Ki p _ u(K, I - 1) - v(K - 1, I) + Dv(K - 1, I) - Ev(K, I - 1) + G[x(K, I) -x(K, I - 1)] - F[x(K, I) - x(K - 1, I)] 

D - I 


( 20 ) 


u(K, I) = u(K, I - 1) + E[v(K, I) - v(K, !-!)] + G[x(K, I) - x(K, I - 1)] 


where 
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D = 0. 5[D(K, I) + D(K - 1, I)] 
E = 0. 5[E(K, I) + E(K, I - 1)] 



F = 0. 5[F(K, I) + F(K - 1, I)] 
G = 0. 5[G(K, I) + G(K, I - 1)] 


The iterative process begins with an estimate of the flow properties (i. e. , u, v, and c 
at the field point x(K, I), y(K, I)) which determine the coefficients A, B, D, E, and F. 
The field point and the flow properties are then calculated from equations (19) and (20). 
New values of the quantities A, B, D, E, and F from which the field point and its prop- 
erties can be recalculated are thus obtained. The iteration continues until a specified 
convergence between successive calculations is reached. 


Basic Characteristic Boundary Point 

Calculation of the basic boundary point (which occurs in the fixed -boundary problem) 
is performed in an analogous manner to the basic field point except for the replacement 
of one of the characteristic equations by boundary conditions. The condition of a solid 
boundary requires that the normal component of velocity at the point x(K, I), y(K, I) 

(fig. 4) be identically zero. 

The boundary curve is specified as a polynomial of third order, that is, 

y w = Rq + RjX + Rgx 2 + R 3 x 3 


Hence, the boundary point x(K, I), y(K, I) is determined by the simultaneous solu- 
tion of the C_ characteristic equation with the equation of the boundary curve. Thus, 


R 3 x(K, I) 3 + R^K, I) 2 + (Rj - B)x(K, I) + Rq - y(K, I - 1) + Bx(K, I - 1) = <fi 


y w (K, I) = Rq + RjX(K, I) + R^K, I) 2 + R 3 x(K, I) 3 


J 


( 21 ) 


The velocity components are obtained from the compatibility condition along the C_ char- 
acteristic and the condition of no flow normal to the boundary. This gives 


u(K, I) = fe 1 - 1) - Ev(K, 1-1 ] + G[x(K, I) - x(K, 1 - 1)] 

1.0 - eE 


( 22 ) 


v(K, I) = eu(K, I) 
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where 


e = 


— = ^ + 2.0 R^K, I) + 3. 0 
dx 


R 3 x(K, I ) 2 


The iteration for the basic characteristic boundary point proceeds in a like manner as 
iteration for the basic field point. 


Shock Field Point Calculation 

Shock field point calculations are inherently different from field point calculations 
in that the conditions known are just upstream of the shock wave. Thus, the oblique 
shock relations must be included in the calculations of the flow field behind the shock 
wave. The shock point x(K, I), y(K, I) is geometrically located on the downstream side 
of the intersection of the shock wave and the opposite family of characteristics in the 
upstream region (fig. 5). The upstream properties of the shock point x(K, I), y(K, I) 
are found by a linear interpolation between the proper net points in the upstream region. 
The iterative process to determine the downstream conditions begins with the assumption 
that the deflection angle at the point x(K, I) , y(K, I) is the same as the deflection angle 
at the previously calculated shock point x(K - 1, I), y(K - 1, I) (fig. 6). From the 
oblique shock relations, the first estimate of the flow properties behind the shock are 
established. A C + characteristic is then constructed from the shock point x(K, I), 
y(K, I), as shown by the dashed line in figure 6. The intersection x re j, y re ^ of this 
characteristic with its opposite member is the reference point. Conditions at this refer- 
ence point are obtained by a linear interpolation between the previously calculated points 
x(K - 1, I), y(K -1,1) and x(K - 1, I + 1), y(K - 1, I + 1). A basic characteristic field 
point calculation, previously described, is then performed to obtain the lattice point 
x(K, I + 1), y(K, I + 1). This calculation is based on the two net points x(K, I), y(K, I) 
and x(K - 1, 1+ 1), y(K - 1, I + 1). The shock point x(K, I), y(K, I) is then recomputed 
by a field point calculation based on the lattice points x re ^, y re j and x(K, I + 1), 
y(K, I + 1). A new estimate of the deflection angle at the shock point is then obtained, 
and the sequence is repeated until the desired convergence has been obtained. 


Shock Boundary Point Calculation 

To obtain the quantities at a shock point x(K, I), y(K, I) under the conditions where 
the previous shock point lies on a physical boundary (fig. 7), a C + characteristic is 
constructed from x(K, I), y(K, I) (dashed line) which intersects the boundary. The 


12 



intersection point of this characteristic with the boundary becomes the reference point 
x ref’ ^ref' T ^ e con ditions behind the shock at the point x(K, I), y(K, I) are obtained in 
the manner described in the previous section. The iteration proceeds in a similar man- 
ner as the shock field point calculation except the properties at the points x re £, y re ^ 
and x(K, I + 1), y(K, I + 1) are obtained by a basic boundary point calculation previously 
described. 


CONSTRUCTION OF FLOW FIELD 

Two basic methods of specifying boundary data were used in the computer program. 
These are (1) the fixed-boundary problem in which the centerbody and cowl contours are 
the primary input variables to the program and (2) the free -boundary problem where a 
single surface curve and Mach number distribution are the primary input data. In the 
actual construction of the flow field, both methods are used in conjunction with each other. 


Cowl Lip Calculation 

One example of a fixed -boundary problem is illustrated in figure 8. The two bound- 
ary curves, y = y c b(x) an< * y = y c (x), represent the centerbody and cowl surfaces of the 
diffuser. The initial data lie along the incident shock; however, these data are incre- 
mentally determined as the flow field behind the incident shock is constructed. The con- 
struction begins with a shock boundary point calculation, the boundary point being located 
at point A in figure 8. The calculations always start with a downstream characteristic 
from the initial datum line, in this case the incident shock. Since the boundary has been 
reached in the first calculation (point A) , the computations proceed with a shock field 
point calculation to locate the next shock point. Field point calculations are then per- 
formed to determine the next downstream characteristic from the shock to the cowl 
boundary, until the point just before the boundary has been reached. At that time, a 
boundary point calculation is performed to determine the flow properties on the cowl 
surface. This sequence of calculations is continued until the entire incident shock has 
been constructed. 

At this point in the calculations, there are several alternatives available in the com- 
puter program. These alternatives are discussed in the following sections. 
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Reflecting Shock Problem 


\he incident shock in figure 8 can be allowed to reflect off the centerbody surface 
represented by y = y c k(x). In this case, the construction of the flow field behind the re- 
flected shock proceeds in exactly the same manner as the flow field behind the incident 
shock, which has previously been discussed. 


Canceled-Shock Problem 

Another variation of the fixed -boundary problem which is used in the computer pro- 
gram is to specify a lower surface y = y c k(x) such that the shock is canceled at that 
boundary. The initial datum line now becomes the last downstream characteristic from 
the incident shock to the upper surface, y = y (x). Under these conditions, the lattice 
points along the initial datum line are known. Hence, only field point and boundary point 
calculations are used sequentially in the construction of the flow field. If two character- 
istics of the same family intersect, the upstream characteristic is deleted and the cal- 
culations continue. Under this condition, the program automatically prints out the locus 
of the intersecting characteristics of the same family. 


Free-Boundary Problem 

Another alternative which is provided in the computer program is to specify that the 
incident shock (fig. 8) be canceled and to furnish the centerbody contour, y = Y c ^,(x), 
and the velocity distribution along this surface. This is the free-boundary problem and 
is illustrated in figure 9. In this case, the initial datum line is specified to be the first 
downstream characteristic from the upper surface that intersects the lower surface of 
the diffuser (line AB in fig. 9) which was determined from the previous calculations. 

The calculations begin with a stream function integration along the initial character- 
istic AB from the lower to the upper surface (see appendix C). This integration estab- 
lishes the numerical value of the streamline which is the upper surface \j/ = ^ ref since 
the point A lies along this contour. The numerical value of the streamline which makes 
up the lower surface is arbitrary and, hence, set equal to 1 (i// = 1). The construction 
of the flow field and the determination of the unknown surface, represented by ip = 4* re f> 
proceed in the following manner. Since the contour and velocity distribution along the 
lower surface are specified, the sequence of calculations begins with a field point cal- 
culation for the second lattice point D along the characteristic CE. Simultaneous 
with the field point calculation, a mass flux integration is performed. This se- 
quence of calculations is continued along the characteristic CE until the calculated 
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stream function first exceeds the value of the reference stream function (i. e. , \p = ^ re f)» 
At this time, an iteration is performed to locate the point along the characteristic CE 
for which the value of the stream function matches the reference value (i. e. , \p = ^ re f) • 
This iteration is based on the two lattice points < ty re f and i// > ip re f indicated by 
square symbols in figure 9. The boundary coordinates (point D) in addition to the prop- 
erties at this point are thus established by this iteration. The calculations are continued 
in the manner described until the entire flow field has been established. 

Another example where specifying the contour and the velocity distribution proves 
advantageous is in the construction of the flow field resulting from an isentropic spike 
(fig. 10). In this example, the upper surface degenerates to the focal point; hence, this 
point is a singularity. Conditions at the focal point are specified by assuming two- 
dimensional reverse Prandtl -Meyer flow. The initial data are obtained from conical 
flow field calculations at specified intervals. The reference stream function ip re f is 
specified by integration of the mass flux along the initial characteristic from the focal 
point to point A on the cone surface. Field point calculations are then performed simul- 
taneously with a stream -function integration to establish the next leading characteristic 
through the focal point, OB in figure 10. When the value of the stream function first 
becomes less than the value of the reference stream function t// < an iteration is 
performed to locate the surface coordinates, point B. Successive characteristics are 
then constructed until the flow field and spike contour have been established. 

In principle, the establishment of one physical boundary by means of a mass flux 
integration does not preclude flows where shock waves are present. Under these condi- 
tions, however, the integration must account for the change in entropy experienced 
across the shock wave. 


SELECTED INLET EXAMPLES 

To indicate the use of the computer program in the design of supersonic inlets, four 
examples are presented and discussed. The purpose of this section is to indicate the 
manner in which the computer program can be used to satisfy different design require- 
ments. There are two types of calculations which can be used to solve for the internal 
flow field. In one type the inlet geometry is specified (fixed -boundary problem) and in 
the other the Mach number and surface angle distribution must be specified (free- 
boundary problem) . In addition, additive drag and mass flow spillage calculations are 
included (see appendixes D and E) for conical flow or simple wedge spillage. For those 
inlet configurations for which the boundary is to be determined, the computer program 
will curve-fit the unspecified surfaces, and the coefficients will appear in the output 
listings. 
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This section also illustrates the use of the primary subroutines and the functional 
organization of the main program. FORTRAN listings of the main programs and sub- 
routines are found in appendix F. A complete listing of the program is presented in 
appendix G. 


Bicone Inlet 

Presented in figure 11 is a 10°-18. 5° bicone mixed -compression inlet designed for a 
free-stream Mach number of 2. 50. The second oblique shock intersected the cowl lip 
for the design condition (fig. 11(b)) while about 0. 5 percent of the capture mass flow was 
spilled supersonically by positioning the cone shock forward of the cowl lip. The initial 
internal cowl angle was set at 5. 0°, while a Mach number of 1. 30 and nominal flow angle 
of -1. 4° were specified for the throat conditions. In addition, it was specified that a 
constant throat Mach number section extending a distance downstream of the throat sta- 
tion be included in the inlet design. Under these conditions, about 60 percent of the 
supersonic area contraction took place externally. The theoretical total -pressure re- 
covery at the throat was 0. 968 behind the terminal shock. 

The conical shock angle and boundary conditions on the first cone were computed in 
subroutine CSA by specifying a free-stream Mach number (AMO) of 2. 500 and a cone 
angle (THETAC) of 10. 0°. The profile of the second shock as well as the resulting flow 
field were computed in subroutine CONE. This was accomplished by specifying the de- 
flection angle (ALPHA) between the two cones and the initial spacing parameter (START) 
which locates the original net -point spacings along the shock. The coordinate system 
always places the origin at the spike tip. For the example shown in figure 11, the param- 
eter (START) was set to a value of 0. 050, which gave a total of 27 net points along the 
second shock. The cowl lip is located at unit radius by specifying the angle between the 
centerline and angular cowl position (THETAL) and setting SPILL = 1.0. This positions 
the intersection of the two cones relative to the cowl lip by requiring the second oblique 
shock to intersect the cowl lip. Additive drag and mass flow spillage calculations are 
then performed in subroutine DRAG, assuming only conical flow spillage (see appendix E). 

The internal flow is constructed by using subroutines SHOCK and SURF in combina- 
tion. Subroutine SHOCK constructs the cowl oblique shock and the resulting flow field by 
using the specified initial cowl angle (COWLA), which in this example is 5.0°. If there 
are no coefficients specified for the cowl surface, the program automatically assumes a 
straight line whose slope is tan(COWLA) and which passes through the cowl lip position. 
For this particular example, the cowl oblique shock was canceled at the centerbody sur- 
face, x = 2. 78 in figure 11(b). This was accomplished in the program by setting 
NSHK = 1 (i. e. , there should be only one oblique shock internal to the inlet with no shock 
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reflections). Under this condition, subroutine SHOCK locates the first characteristic 
behind the cowl oblique shock which intersects both the cowl and centerbody surfaces and 
performs a mass flux integration to determine the reference streamline. This estab- 
lishes the initial datum line to be used in subroutine SURF and the reference streamline 
of the cowl. All other information computed downstream of this initial characteristic is 
thus ignored. Hence, under the condition where NSHK = 1, subroutine SHOCK con- 
structs the flow field within the domain in which the cowl surface can influence the 
oblique shock profile. Downstream of the initial characteristic which intersects the 
centerbody determined in subroutine SHOCK, the internal flow as well as the contour 
coordinate points are constructed in subroutine SURF. This is accomplished by specify- 
ing in the main program as input data the parameters AMT(J), THR(J), ANG(J), and 
NIS(J) (see appendix F). With this information, subroutine SURF constructs a parabolic 
surface of length THR(J), measured from the last known centerbody surface point, and 
having a final angle of ANG(J). On this surface, the program assumes a linear Mach 
number distribution which has a final Mach number of AMT(J) at the end of the section. 
Subroutine SURF, with this information, constructs the cowl contour which satisfies the 
boundary conditions specified. The number of net points for this calculation is controlled 
by the index parameter NIS(J). It can be seen that very complex surfaces which have 
precise distributions can be constructed from the basic forms described by specifying as 
many sections (indicated by the J index) as needed. There are several options con- 
structed into this program which serve useful functions. By setting THR(J) = 0, the 
program will focus all leading characteristics from the cowl at the point which starts 
the J section. Conditions at the focal point are automatically assumed to satisfy the 
Prandtl-Meyer relation. By setting ANG(J) = 0, subroutine SURF will construct a 
centerbody surface whose Mach number distribution is specified by AMT(J) and whose 
angular distribution is satisfied by focused Prandtl-Meyer compression at a point below 
the body surface. The final angle in this case is determined by the program. For the 
example shown in figure 11, the input data were as follows: 


SEGMENT 

AMT(J) 

THR(J) 

ANG(J) 

NIS(J) 

1 

1. 300 

0.400 

0 

10 

2 

1. 300 

.400 

0 

8 

3 

1.300 

.400 

0 

8 


Subroutine SURF, in this example, constructs the internal contours of the inlet such that 
the flow is initially compressed isentropically to a Mach number of 1. 300 over a dis- 
tance equal to 40 percent of the cowl lip radius, measured along the centerline from the 
shock impingement point. This segment of the centerbody surface was assumed to be a 
streamline constructed from the Prandtl-Meyer relations. Since the flow on the center- 
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body has already been compressed to Mach 1. 300, the data input for segment 2 was used 
to construct a transition section on the cowl to reduce its Mach number to 1. 300. This 
transition section on the cowl can be seen in figure 11 and lies between x = 3. 02 and 
x = 3. 22. The data for segment 3 were used to hold the 1. 300 throat Mach number a dis- 
tance of 0.400 downstream of segment 2. When the inlet flow field calculations are com- 
pleted, the program curve-fits the centerbody and cowl contours and lists the table of 
coefficients as output. These coefficients can then be used to evaluate the off -design 
capability of the inlet. 


Two-Dimensional Focused Compression Inlet 

Shown in figure 12 is a two-dimensional isentropic ramp inlet, designed for a free- 
stream Mach number of 2. 70, in which the external flow was compressed to a focal point 
Mach number of 2.05. The initial ramp angle and internal cowl angle were both set at 
5. 0°, while at the inlet throat a Mach number of 1. 300 and flow angle of -5. 00° were 
specified. With these conditions, approximately 70 percent of the supersonic area con- 
traction took place externally. The Mach number distribution and the characteristic 
solution along the centerbody and cowl contours are presented in figures 12(a) and (b), 
respectively. For this inlet, the theoretical total -pressure recovery behind the terminal 
shock was computed to be 0. 956. 

For this inlet example, it was specified that two-dimensional flow calculations 
should be performed. This was accomplished by setting NDIM = 2, as opposed to 
NDIM = 3 for the previous example. The calculations to construct the isentropic ramp 
contours and flow field were performed in subroutine SPIKE. Again, the spacing param- 
eter START must be specified to establish the net -point spacings along the initial datum 
line. In addition, an M-index parameter must be specified to divide the flow field into 
M + 1 C + characteristics passing through the focal point. Thus, the M-index parameter 
divides the Mach number distribution at the focal point into M equal Mach number seg- 
ments between the initial Mach number and the focal point Mach number (AMF). Down- 
stream of the isentropic section, the flow field is constructed in subroutine SHAPE. The 
initial datum line in this subroutine is taken to be the last of the family of C + character- 
istics computed in subroutine SPIKE. Unless otherwise specified, the program auto- 
matically assumes a straight line segment aft of the isentropic ramp, as was the case 
for the example shown in figure 12. This condition is indicated by the constant Mach 
number distribution along the ramp surface downstream of x = 1. 34 (circular symbols 
in fig. 12(b)). 

For this inlet example, the cowl oblique shock was again canceled at the ramp sur- 
face. The internal flow field and surface contours were determined in the same manner 
as for the bicone inlet previously discussed. 
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Two-Dimensional Ramp Inlet 


Another example of a two-dimensional inlet that was designed by the computer pro- 
gram is shown in figure 13. Although similar in appearance to the example previously 
discussed, basically different design specifications required a somewhat different com- 
putational approach. Foremost in the design specifications was the condition that the 
initial cowl angle be set at 0°. For optimum overall performance, this requirement 
dictated that, at most, 50 percent of the supersonic area contraction takes place exter- 
nally. The throat Mach number and flow angle were set at 1. 300 and -5.00°, respec- 
tively; the same as the previously discussed two-dimensional inlet. To accomplish 
these objectives, the external ramp contour was designed such that all the compression 
waves entered the inlet (fig. 13(b)). This had the desirable effect of reducing the de- 
flection angle at the cowl lip from the 12. 2° of the previous example to 5. 0°, while keep- 
ing the same amount of flow compression on the ramp surface upstream of the oblique 
shock. Consequently, the theoretical total -pressure recovery behind the terminal shock 
in the throat was 0. 971 for a free-stream Mach number of 2. 700. 

The external contour of the ramp surface was constructed in two segments. The 
first segment was a plane two-dimensional ramp with an angle of 5. 0°. The second seg- 
ment was composed of a parabolic section chosen such that the final Mach number on the 
ramp surface before the cowl oblique shock was about 2. 00. Computational work for the 
external flow field was performed in subroutine CONE by specifying ALPHA = 0 and 
introducing the proper coefficients for the ramp surface. The internal contours of this 
inlet were determined in the manner described in the previously discussed examples. 


Self-Starting Inlet 

The condition of self -starting with a minimum cowl wave drag imposes unusual com- 
putational requirements on the fourth inlet example (fig. 14) so that it warrants discus- 
sion. Externally, the flow was compressed isentropically to a final focal point Mach 
number (AMF) of 1. 600. The flow field and contour of the isentropic spike were com- 
puted in subroutine SPIKE. The initial half -cone angle (THETAC) was set at 16. 0° and 
the free-stream Mach number was specified to be 2. 500. It was also specified that for 
an average entering Mach number of 1. 60, the minimum internal Mach number should 
not be less than 1. 35. In addition, the requirement of minimum cowl wave drag dictated 
that the initial internal cowl angle be about 17.0°, which set the Mach number just behind 
the oblique shock on the cowl at 1. 35. In order to turn the cowl contour back as quickly 
as possible to reduce the cowl projected area, external expansion was provided on the 
spike contour aft of the station x = 1. 61. The contour of this segment of the spike was 
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chosen to provide a nearly uniform Mach number profile at the cowl lip. Results for 
this section of the flow field were computed in subroutine SHAPE by specifying the ap- 
propriate coefficients. A nominal throat flow angle of -3. 0° was set for this inlet con- 
figuration. The problem, therefore, in this inlet design reduces itself to finding the 
proper Mach number distribution and contour along the centerbody surface which would 
ensure a constant Mach number of 1. 35 along the cowl surface while turning the flow 
along the cowl a total of 20. 0°. 

The calculation for the solution of this problem was accomplished in subroutine SURF 
by systematic variations of the program input variables AMT(J), THR(J), and ANG(J). By 
successive approximations, several inlet designs were established, each of which came 
nearer to satisfying all the requirements. The characteristic solution for the final con- 
figuration is shown in figure 14(b), while the Mach number distribution is shown in fig- 
ure 14(a). The overall total -pressure recovery behind the terminal shock in the throat 
region was 0. 955 for a design Mach number of 2. 500. 


Off-Design Calculation 

Off -design calculations can be performed on the various types of inlet configurations 
discussed. In each case, the coefficients of the equations specifying the various surfaces 
are input data to the program. If the cowl oblique shock was to be canceled at the center- 
body shoulder for on-design operation, as was the case in the inlets discussed, the loca- 
tion of the shoulder must be specified for off -design calculations. This is done by 
specifying the index parameter (ISHK) which selects the position point BODY(ISHK) as 
the shoulder point. It must be understood that canceling the shock for design conditions 
does not necessarily imply cancellation at off -design Mach numbers. 

The position of the cowl lip for off -design calculations is specified by the coordinates 
X COWL, Y COWL. If these two parameters are nonzero quantities, the program will 
automatically shift the coefficients describing the cowl contours to their proper positions. 
If the cowl lip is positioned outside the region of influence of the previously computed 
flow field, the calculations will terminate, and an error message will be indicated in 
the output listings. This can be changed by choosing another value of BODY(2) along 
the centerbody surface which initiates the forward flow field calculations. 


COMPARISON WITH EXPERIMENTAL DATA 


To illustrate the agreement between the numerical methods described in the previous 
sections and experimental data, selected examples are presented in figures 15 to 18. 
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Attention is given in the discussion of these examples to show the influence of boundary 
layers and to indicate the degree of accuracy of this type of flow calculation. 


Isentropic Spike Configurations 

Figure 15 presents an isentropic spike forebody configuration designed for a free- 
stream Mach number of 2.49 and a final focal point Mach number of 1. 46. The char- 
acteristic solution is shown in figure 15(b), while a comparison between the measured 
and calculated static -pressure distribution is presented in figure 15(a). The previously 
unpublished data were obtained from the investigation of reference 6. In general, the 
boundary layer did not greatly influence the static -pressure distribution along the spike 
surface, as testified by the good agreement between the calculated (solid line) and 
measured (symbols) values. Some boundary-layer influence was noticed, however, near 
the end of the isentropic section, but the focusing characteristics were not appreciably 
impaired, as indicated by the Schlieren photograph (fig. 15(c)). 

More pronounced boundary-layer influence was indicated on an isentropic spike 
forebody configuration designed for a free-stream Mach number of 3. 85 (fig. 16). The 
axisymmetric characteristic solution for this configuration is shown in figure 16(b). The 
cone had an initial half angle of 8. 0°, while the final focal point Mach number was set at 
1. 75. Figure 16(a) shows a comparison between the calculated (solid line) and measured 
(symbols) static -pressure distribution along the spike contour obtained from reference 7. 
Generally, good agreement was obtained between the experimental and calculated static 
pressures; however, all data tended to fall somewhat higher than the predicted distri- 
bution. This was typical everywhere except for a short range in the no-roughness condi- 
tion (circular symbols) immediately following the laminar boundary -layer separation 
region (dashed line). Turbulent boundary layer was undoubtedly induced along the spike 
surface by the application of roughness at the spike tip, as indicated by the slightly 
higher static -pressure level (comparing circular and square symbols). This was 
realized everywhere except in the separation zone of the no-roughness case. 


Bicone Forebody Configuration 

Shown in figure 17 is a bicone forebody configuration in which the first cone had a 
half angle of 10. 0° while the second cone half angle was 18. 5°. Calculations were per- 
formed at a free-stream Mach number of 2.49, which was identical to that for the data. 
The data were obtained in an unpublished investigation conducted in the Lewis 10- by 
10 -foot supersonic wind tunnel at maximum Reynolds number. The diameter of the 
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oblique shock intersection point was 40. 6 centimeters. The axisymmetric characteristic 
solution is presented in figure 17(b). Coordinate points along the contour were normal- 
ized such that the intersection of the first and second oblique shocks occurred at a radius 
(dimensionless y coordinate) of 1. A comparison between the experimental data (cir- 
cular symbols) and calculated values (solid line) of static pressure along the surface con- 
tour indicated generally fair agreement (fig. 17(a)). Upstream of the junction between 
the first and second cones, the data were slightly lower than calculated values. Down- 
stream of the juncture point, slightly higher static pressures were measured than were 
calculated, indicating some boundary -layer influence. Surface disturbances on the 
second cone, clearly visible in the Schlieren photograph (fig. 17(c)), prevented the pos- 
sibility of distinguishing any boundary-layer influence aft of x = 1. 40. There appears to 
be no boundary -layer "bridging" at the junction of the conical surfaces. This type of 
flow separation may induce a greatly modified shock pattern, particularly in the vicinity 
of the cowl lip of an inlet. 


Single-Cone Inlet Configuration 

Figure 18 presents the computer solution obtained for a single-cone mixed- 
compression inlet designed for a free-stream Mach number of 2. 500. External com- 
pression was obtained by a 12. 5° half -angle conic forebody. The initial cowl angle was 
set at 0° in this inlet example; the oblique shock originating from the cowl lip was re- 
flected from the internal contour surfaces. The internal contour surfaces were obtained 
by curve-fitting the contour specifications used in fabricating the inlet. 

The characteristic solution for this inlet configuration is presented in figure 18(b). 
Shown in figure 18(a) is a comparison between the static -pressure distribution obtained 
from the computer program (solid and dashed lines) and experimental data (symbols) 
obtained from reference 8. In general, there was good agreement between the theoreti- 
cally determined static pressures and measured pressures upstream of the throat re- 
gion (x = 3.48). Within the throat section, agreement between calculations and data is 
more difficult to realize owing to the complex shock boundary -layer interactions. These 
interactions were partially minimized by porous bleed in the regions indicated in fig- 
ure 18(b). In general, the shock reflection points (as indicated by data) appear to occur 
upstream of the calculated points and thus indicate reflections off the boundary layer. 

This tends to foreshorten the shock regions in the downstream direction. More detailed 
discussions of the experimental data can be found in reference 8. 

Off -design performance calculations of this inlet at Mach numbers of 2. 30 and 2.02 
reveal increasingly better agreement in the throat region with decreasing free-stream 
Mach number. This is probably due to the lessening boundary -layer influence. 
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CONCLUDING REMARKS 


The theory, numerical methods, and computer program for the design of supersonic 
inlets have been presented. The computer program was written with the intention of de- 
signing many types of inlets with a minimum of computer time. For example, the execu- 
tion time for the inlet presented in figure 11 was 0. 62 minute with 30 initial net points 
along the starting shock. The execution time for the inlet presented in figure 18 was 
1.23 minutes for five internal shocks with 31 net points along the first cowl oblique shock 
wave. 

Because computer programs of this type are used primarily for design application, 
it becomes important for the programmer to use methods whereby optimum designs are 
more readily obtained. By introducing throat Mach number and flow angle as direct pro- 
gram input variables, the difficulty in obtaining the proper surface contour was consider- 
ably reduced. This reduced the overall time required to obtain the designed inlet con- 
figuration. In addition, minimum inviscid throat Mach number distortions could be 
achieved. 

Lewis Research Center, 

National Aeronautics and Space Administration, 

Cleveland, Ohio, September 3, 1968, 

126-15-02-11-22. 
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APPENDIX A 


DERIVATION OF GENERAL COMPATIBILITY 
AND CHARACTERISTIC RELATIONS 

The set of equations which describe steady, irrotational isentropic flow in three 
dimensions with cylindrical symmetry are 

2 

(c 2 - u 2 )u x - uv(Uy + v x ) + (c 2 - v 2 )v y + — = 0 (Al) 

v x ” Uy = 0 (A2) 

ct 2 (u 2 + v 2 ) + (1 - ct 2 )c 2 = 1 (A3) 

where the velocity components u,v and the speed of sound c are dimensionless with 
respect to the critical speed, and the coordinates of a point x,y are dimensionless with 
respect to a characteristic length. In equation (Al) the subscripts x,y indicate differ- 
entiation. A linear combination of equations (Al) and (A2) is 

pu x + qUy + rv x + sVy + t = 0 (A4) 

where 

p = Xj(c 2 - u 2 ) 
q= -X 2 

r = X 2 - 2XjUV 

s = Xj(c 2 - v 2 ) 

X 1 c 2 v 
t = — 

y 

The characteristic equations are obtained by requiring the derivatives of u,v with 
respect to x,y to form derivatives in the same direction (i. e. , dx:dy = p:q = r:s. If 
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x(t) and y(r) represent a curve with x z :y z = p:q = r:s, equation (A4) represents the 
derivatives of u,v along this curve, and r(x,y) define the characteristics. Thus, 


or 


q . s y T 

py T - qx T = o 

ry T - sx T = 0 


(A6) 


Using the set of equations (A6) gives the pair of linear algebraic equations 

(c 2 - u 2 )y T Xj + x t A 2 = 0 
[-2uvy T - (c 2 - v 2 )x t ] Xj + y T A 2 = 0 

For the set of equations (A7) to have a nontrivial solution for X, the determinant of the 
coefficients of X must vanish; thus, 

(c 2 - u 2 )y 2 + 2uvx T y T + (c 2 - v 2 )x 2 = 0 (A8) 

The characteristic directions are therefore given by 



Cx = — = - 


X T C 2 -U 2 


-uv ± V(uv ) 2 - (c 2 - u 2 )(c 2 - V 2 ) 




(A9) 


and the characteristics equations become 


y a = C + x a 

^ = ^- x /3 


(A10) 


Equations (A10) are two separate ordinary differential equations of first order which 
define two one-parameter families of characteristics C + and C_ in the x,y plane. 
These families are represented in the form a(x,y) = Constant, /3(x, y) = Constant and 
form a curvilinear coordinate net such that j3 is constant along a C + characteristic 
and a is constant along C_. To obtain the compatibility conditions, equation (A4) was 
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successively multiplied by x and y. Using equation (A6) gives 


pu^ + nr + tx„ = 0 

T T T 


qu + sv^ + ty^ = 0 

n T T J T 

Using the set of equations (A5) results in the pair of algebraic equations 




(c 2 - u 2 )u t - 2uvt t + — X T 


h + V2 = 0 


(c 2 -v 2 )v T + ^y T 

y 




X 1 " U t X 2 = 0 


(All) 


(A12) 


Using the first equation of sets (A7) and (A12), respectively, and requiring that this pair 
has a nontrivial solution for X gives the results 


u T - c± v T - ♦ l-i— xV. = 0 


: 2 -u 2 / T \c 2 - u 2 V 7 


The pair of compatibility conditions are thus given by 


~\ 


u a + ?_ v a + ■ — = 0 


2 2 y 

i - u y / 


"a 


+ ?+Vp + 1^7 2 y r °. 


(A13) 


The characteristic equations and compatibility relations can be more conveniently ex- 
pressed by introducing the angle S between the flow direction and the positive x-axis; 
then, 


u = q cos 0 
v = q sin S 


(A14) 


The roots ? + and £_, being the slopes of characteristic directions of C + and C_, are 
thus 
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K + = tan (a + m) 
£ _ = tan (a - m) 


(A15) 


where M is the angle between the flow direction and that of the corresponding character- 
istic. 

The characteristic expression, equation (A8), can be written in the form 

c 2 (l + e ± 2 ) - (u? ± - v) 2 = 0 
Using equations (A14) and (A15) gives the following result: 

c 2 = q 2 sin 2 p (A 16) 


Thus, the speed of sound is simply the component of flow velocity normal to the direction 
of a characteristic. 

With these rotations, the characteristic equations assume the form 


Y a = tan(a + p)x a C + 
y^ = tan(a - p)x C_ 


u = -tan(a - p)v - — x C 

01 2 2 + 


c - u 




U£ = -tan(a + p)v /3 


2 2 
c - u 


x C 




(A17) 
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APPENDIX B 


DERIVATION OF DIFFERENTIAL EQUATIONS OF CONICAL FLOW 

For the mathematical construction of conical flow patterns, assume that the velocity 
components u,v and hence the speed of sound c depend only on the ratio 


x 

V = — 

y 


(Bl) 


Equation (A2) thus becomes 


v + mi =0 
V 1 V 


(B2) 


and equation (Al) becomes 


(c 2 - u 2 )u - 2uvv - (c 2 - v 2 )r 7 V + c 2 v = 0 
' ' r\ V V 


(B3) 


where subscripts refer to differentiation. Equation (B3) assumes a particularly useful 
form when v is introduced as a function of u; thus, from (B2) 


n 


= -JL=-y 


u 


u 


(B4) 


Differentiating (B4) with respect to u leads to the following relation: 


u v 


uu 


V = 

V 


u 


uu 


(B5) 


Introducing equations (B4) and (B5) into equation (B3) leads to the particularly simple 
form of the Taylor-Maccoll equation 


, 2 - 2 / \2 
vv uu = 1 + v u - c (u + vv u } 


(B6) 
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Eliminating c 


2 


by mean of Bernoulli's equations yields 


vv. 


uu 



(1 - cr 2 ) (u + w u ) 2 

i 2/2 2 X 

1 - <t (u + v ) 


(B7) 


Along the cone surface, the flow has the direction of the ray rj = x/y traced by the cone 
in the x,y plane, and thus 


V = 


u 




> 


V = - — 

u V 


(B8) 


The conditions to be satisfied along the conical shock (ref. 1) are given by 


u = (1 - or 2 )qQ cos 2 /3 g + 

q 0 

v = (q 0 - u) cot /3 g 


(B9) 


where is the dimensionless free-stream velocity and j3 g is the conical shock angle. 
In addition, the initial slope along the shock is given by 


v u = 


u - 


% 


(BIO) 
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APPENDIX C 


DERIVATION OF STREAMLINES IN AXISYMMETRIC FLOWS 


The continuity equation for axisymmetric flow may be written as 


(puy) x + (pvy) y = 0 


From equation (Cl), there exists a stream function x,y such that 


= _pvy 


= puy 


(Cl) 


(C2) 


Along the two one -parameter families of characteristics C + and C_ in the x,y plane, 
which are represented by j3(x, y) = Constant and a(x, y) = Constant, 


^x x a + ^y^a 

*/3 = V/3 + >V/3 
where the derivatives x a , y a and x^, y^ are given by 

x a = cos (5 + ix) along C+ 




y a = sin (3 + n) along C + 

Xq = cos(d - /it) along C _ 

y^ = sin(3 - ju) along C _ 

When cylindrical coordinates are introduced such that 

u = q cos 3 
v = q sin S 


> 


(C3) 


(C4) 


and the relations, equations (C4) are used, equation (C3) becomes 
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>p a = qy [sin (3 + p) cos 3 - cos(3 + p) sin 3 ] ''j 
^ - qy [sin(3 - p) cos 3 - cos(3 - p) sin a] j" 


(C5) 


Simplifying the above equations results in 


i K, = qy sin M = 

M 


"\ 


^ = -qy sin M = 


_ -pqy 

M 


From the set of equations (C4), 


do = 


dy 


sin (p + 3) 

d/3=-_ iz 

sin(ju - 3) 


(C6) 


(C7) 


Therefore, 


*B = ^A + 



pqy dy 
M sin(p ± 3) 


(C8) 


where the + sign is used for integration along a C + characteristic and the - sign is 
used for integration along a C_ characteristic. If the interval AB is used as the dis- 
tance between successive points on the characteristic net (along the appropriate char- 
acteristic), it maybe assumed that pq/M sin(p ± 3) takes the average of the values at 
points A and B; thus, 



For two-dimensional flows, equation (C9) reduces to the expression 

= + f — "I + r — 1 i (y B " y A* ( cl °) 

2 |1_M sin(p ± 3)J a [M sin(p ± 3)J B J 


31 


The derivation of the stream -function integration procedure along a C + character- 
istic appeared first in reference 9. For a simple wave region (two-dimensional flow) 
the integration of equation (C 8) takes on a particularly simple form since the integrand 
becomes constant along either the C + or C_ family of characteristics; thus, 


*B = *A + 


pq 

M sin(/Lt ± 5) 


(y B - y A > 


(Cll) 
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APPENDIX D 


DERIVATION OF STREAMLINES IN CONICAL FLOW FIELD 

The construction of streamlines in a conical flow field becomes particularly simple 
if the integration of mass flow is performed along rays passing through the apex of the 
cone. Thus, equation (C3) is replaced by 

° V r + Vr 

where 

r = Vx 2 + y 2 

The derivatives x r and y r are given by 

x r = cos © 
y r = sin © J 

where 0 is the ray angle. Using equations (C2) and the polar form of the velocity com 
ponents gives 


(Dl) 


(D2) 


i// = pq y(sin 0 cos 5 - cos 0 sin 3) 


= pq y sin(© - 3) 


(D3) 


Therefore, when equation (D2) is used, the streamline is specified by 

/•y b 

* = * A + I j ffl y dy <D4) 

J y A 

Since integration of the mass flux is performed along rays passing through the cone apex, 
the quantity inside the brackets is constant, and hence 
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When 
of y 


I 


^ + sin ~ - - - ( y l ■ y i) 

(_ sin 0 J 2 ' ' 


^A» y A = 0 with no loss in generality, the stream function becomes a function 
along a specified ray of angle 0; thus, 


<// = P(Q)q( Q ) sin[0 - 3(0)] y 2 
2 sin © 


(D5) 
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APPENDIX E 


ADDITIVE DRAG AND MASS FLOW SPILLAGE IN CONICAL FLOW FIELD 

For a conical flow field, the transition across the shock wave is governed by the 
oblique shock relations and is followed by a continuous isentropic compression to surface 
conditions (fig. 19). The capture streamline is established by an integration of the mass 
flux along successive generatrix rays lying between the cowl lip location (point P in 
fig. 19) and the shock wave. Derivations of the integral formulas used in this calcula- 
tion are presented in appendix D. The reference stream -function value (capture stream- 
line) is first established by integration of the mass flux along the generatrix ray OP by 
using equation (D5); that is, 




ref 


_ P(O p )q(Op) sin 0 p - 3(© p ) 


2. 0 sin 0„ 


y(© p )‘ 


(El) 


Subsequent points along the capture streamline (i. e. , points Q, R, and S in fig. 19) 
are thus obtained from the relation 


y(©) = 


i 


2.° ^ ref (sin 0) 
p(0)q(0) sin 0 - $(©) 


(E2) 


When the generatrix ray becomes coincident with the conical shock wave, the capture 
mass flow is obtained from the equation 


m _ W (0 S )2 _ y(0 S )2 

m ° Po q o y(0 p) 2 y( 0 p ) 2 

When the notation in figure 19 is used, the additive drag coefficient is defined as 


r Y s 

C d,a J o I (P-P0 )ydy 

rp o Mjy(0 p )^ •'Yp 


(E4) 


Equation (E4) may be approximated by the method of numerical integration. Thus, the 
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additive drag was computed by using the relation 


C 


d,a 


1 

yM§y(0 p ) 2 



P(j j 1) 
P 0 


- 2 


[y(j) 2 - y(j + i) 2 ] 


(E5) 


where the J-index designates the intersection points of the capture streamline and gen- 
eratrix rays, starting from the cowl lip position. 
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APPENDIX F 


COMPUTER PROGRAM FOR DESIGN OF SUPERSONIC INLETS 

The required input data are presented for this program so that it may be readily em- 
ployed. In actuality, there are two main programs, INLET 1 and INLET 3, to assemble 
the subroutines to arrive at the desired inlet design. The difference between INLET 1 
and INLET 3 is the type of forebody configuration desired for the inlet. INLET 1 designs 
an inlet with a focused isentropic forebody configuration, while INLET 3 assumes a bi- 
cone or biramp forebody. The input data required for both of these main programs are 
the same except for three parameters which are noted. In addition, the computer pro- 
gram contains several options, depending on the type and amount of calculations required. 
This is made clear in the following outline of the input required for this program. 

INPUT VARIABLES AND EXPLANATION 


Card 

FORTRAN 

symbol 

Card 

columns 

Description and comments 

1 

AMO 

1-12 

free -stream Mach number 


GAM 

13-24 

ratio of specific heats 


THETC 

25-36 

initial cone or ramp angle, radians 


BETAE 

37-48 

estimate of cone shock angle, radians 


AMF 

49-60 

for INLET 1 only: final focal point Mach number for 
isentropic forebody configuration 


ALPHA 

49-60 

for INLET 3 only: deflection angle 


DELB 

61-72 

initial increment used in calculation for cone shock 
angle; set DELB = 1.0x10"^ 

2 

DELU1 

1-12 

integration increment for conic flow field calculations; 
set DELU1 = 1.0X10 -5 


ERROR 

13-24 

A 

convergence parameter; set ERROR = 1.0x10 for 
most applications 


COWLA 

25-36 

initial cowl angle, radians 


THETAP 

37-48 

for INLET 1 only: location of focal point relative to 
spike tip (measured in radians) 
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Card 

FORTRAN 

symbol 

Card 

columns 

Description and comments 

2 

THETAL 

37-48 

for INLET 3 only: location of cowl lip relative to spike 
tip (measured in radians) 


START 

49-60 

initial net -point spacing parameter 

Both THETAP and THETAL are angular locations relative to the inlet centerline. These 
parameters are used only when the cowl contour is unspecified. Under this condition, 
the program always takes the cowl lip radius as 1. 

3 

FOCUS 

1-12 

for INLET 1 only: set FOCUS = 0 to locate the focal 
point along the initial shock wave; set FOCUS = 1.0 
to locate the focal point at THETAP (In both of these 
cases, the cowl lip is placed at the focal point. ) 


SPILL 

1-12 

for INLET 3 only: set SPILL = 0 to locate the cowl 
along the initial shock wave; set SPILL = 1. 0 to 
locate the cowl at THETAL 


DELP 

13-24 

parameter which limits the extent of the initial datum 
line for off -design calculations (OFF = 1.0); ratio 
of angular location of last initial data point to initial 
shock angle 


The program constructs the initial datum line from the point on the surface, BODY (2), 
to a point in the field whose angular location is DELP times the initial shock angle. 


3 

OFF 

25-36 

set OFF = 0 for on-design calculations; set OFF = 1.0 
to indicate off -design calculations 


XCOWL 

37-48 

x-coordinate of cowl lip position 


YCOWL 

49-60 

y -coordinate of cowl lip position 


The parameters XCOWL and YCOWL are used when a relative shift between the spike tip 
and cowl lip is desired. Under these conditions, it is recommended that the cowl con- 
tours corresponding to the design conditions be used. The program will automatically 
shift the cowl contour such that the cowl lip is located at XCOWL, YCOWL. If there is 
no shift, set XCOWL, YCOWL = 0. If the cowl lip falls outside the region of influence 
of the previously computed flow field, an error message will be printed out and the cal- 
culations terminated. This can be corrected by choosing a different BODY (2) position, 
which is a dummy point used to start the initial flow field calculations. 
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Card 

FORTRAN 

symbol 

Card 

columns 

Description and comments 

3 

PRINT 

61-72 

set PRINT = 0 to bypass the output printout for the 
forebody flow field calculations; otherwise set 
PRINT » 1.0 

4 

NDIM 

1-6 

set NDIM = 3 for axisymmetric flow field calculations; 
set NDIM = 2 for two-dimensional flow field calcula- 
tions 


M 

7-12 

number of initial I-net point spacings for isentropic 
spike or ramp calculations 


NR 

13-18 

number of sets of cowl contour coefficients read in as 
input 


NS 

19-24 

number of sets of body contour coefficients read in as 
input 


NTHR 

25-30 

number of internal isentropic compression sections 


ND 

31-36 

number of increments in additive drag calculation 


ISHK 

37-42 

index which specifies at which centerbody contour point 
BODY OSHK) a discontinuity occurs for cowl shock 
cancellation; used with off -design calculations only, 
i. e. , OFF = 1.0; otherwise set ISHK = 1, indicat- 
ing no shock cancellation 


NSHK 

43-48 

number of internal shock waves to be computed; 
NSHK-1 specifies the number of internal shock re- 
flection points; to cancel cowl lip shock, set 
NSHK = 1 

5 

Ha, 1) 

1-12^ 

cowl contour is specified by a third-order polynomial 
1 of form 


Ra, 2 ) 

13-24 



Ra, 3 ) 

25-36 ^ 

Y = Ra, 1) + Ra, 2)X + Ra, 3)x 2 + Ra, 4)x 3 


4) 

37-42^ 

COWLa) specifies the x-location after which the I-set of 


COWLa) 

43-54 

>■ coefficients is valid; hence, the I-set of coefficients 


1*1, NR 

J 

is valid in the region COWLa) ^ X < COWLa+1) 
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Card FORTRAN 

Card 

Description and comments 

symbol 

columns 


5 S(I,1) 

1-12 ^ 

centerbody is specified by equation 

S(I,2) 

13-24 

> y = S(1, 1) + S(1, 2)X + S(I, 3JX 2 + S(1, 4)X 3 

S(I, 3) 

25-36 


S(I,4) 

37-42 

valid in the region BODY(I) < X < BODY(I+l) 

BODY (I) 

43-54 

BODY(I) performs the same function as the parameter 

1=1, NS 


COWL(I) 

tVi 

The coefficients for the NR and NS set are dummy variables and are set equal to zero, 


while COWL(NR) and BODY (NR) specify the region in which the (NR-1), (NS-1) set of coef- 
ficients are valid. If the cowl and centerbody contours are to be computed, set NR = 1 
and NS = 1 and set the coefficients to zero. If a set of coefficients is provided, for ex- 
ample, for the forebody contour, the program will automatically recompute the internal 
contours for the conditions specified. The following set of cards read in information 
pertinent to computing the internal cowl and centerbody contours, otherwise they are 
unnecessary. 


AMT(J) 

1-12 

final Mach number specified for the end of the J tfl sec 

J=l, NTHR 


tion 

THR(J) 

13-24 

length of J™ section 

ANG(J) 

25-36 

ii. 

final surface angle (radians) for end of J section 

NIS(J) 

37-42 

number of surface net points within THR(J) 
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$ I BFTC INLET3 L I S T , REF , DECK , DEBUG 

c V 

C CHARACTERISTIC PROGRAM FOR TWO RAMP INLET, NDIM=2 

C 

C CHARACTERISTIC PROGRAM FOR TWO CONE INLET, NDIM=3 

C 

C AMO=FREE STREAM MACH NUMBER 

C GAM=R AT 1 0 OF SPECIFIC HEATS 

C THETC= I N I T I AL CONE OR RAMP ANGLE (RAOIANS) 

C ALPHA=SECONO CONE OR RAMP DEFLECTION ANGLE (RAOIANS) 

C COWL A=COWL ANGLE (RADIANS) 

C THETAL=LOCATION OF COWL RELATIVE TO SPIKE TIP (RADIANS) 

C BETAE=EST IMATE OF CONE SHOCK ANGLE (RADIANS) 

C R( I,J)=COEFFICIENTS OF COWL CONTOUR SECTIONS 

C S( I,J)=COEFFICIENTS OF CENTERBODY CONTOUR SECTIONS 

C COWL ( I ) =COWL CONTOUR REGIONS 

C BODY( I )=CENTERBODY CONTOUR REGIONS 

C DELUl=INTEGRATION INCREMENT IN CONE FIELD CALCULATION 

C ERR OR CONVERGENCE PARAMETER, SET ERR0R=0.1-3 FOR MOST APPLICATIONS 

C COWLX=X-COORDI NATE OF COWL LIP 

C COWL Y=Y-COORDI NATE OF COWL LIP 

C 

C M=NUMBER OF INITIAL I -NET SPACE 

C NR=NIJMBER dF COWL CONTOUR SECTIONS 

C NS=NUMBER OF CENTERBODY CONTOUR SECTIONS 

C NTHR-NUMBER OF ISENTROPIC COMPRESSION SECTIONS 

C NSHK=NUMBER OF INTERNAL SHOCKS 

C ND=NUMBER OF INTEGRATION INCREMENTS IN ADDITIVE DRAG CALCULATION 

C I SHK=SHOCK IMPINGMENT POINT ON BODY CONTOUR SECTION 

C 

C COMMAND PARAMETERS 

C 

C SET NDIM=2 FOR TWO DIMENSIONAL FLOW CALCULATIONS 

C SET NO I M = 3 FOR THREE DIMENSIONAL FLOW CALCULATIONS 

C SET AL PHA=0 . 0 FOR A SINGLE CONE INLET 

C SET SP I L L =0 . 0 FOR COWL LOCATED ON FIRST OBLIQUE SHOCK 

C SET SPILL=1 .0 FOR COWL LOCATED AT THETAL 

C SET NSHK=1 FOR ISENTROPIC INTERNAL COMPRESSION 

C SET NSHK GREATER THAN 1, ISHK=I, FOR INTERNAL SHOCKS TO REFLECT 

C SET 0FF=0.0 FOR ON-DESIGN CALCULATIONS 

C SET OFF=1.0 FOR OFF-OESIGN CALCULATIONS 

C 

COMMON AMO, GAM, THETC, DELB, DELU1 ,DELU,UC, VC, DELC,BS, DELS, ERROR, 

1 0, AM, AMU, THETA, A , B , C , D, E , F , G , PRES , DENS , ARE A , DYNP , CP , PS I , 

2 XK,YK ,UK,VK,XJ,YJ,UJ,VJ,XA,YA,UA,VA,DELA,BETAK,SIGK,SIGSO, 

3 X ( 50,50 ) ,Y (50,50) ,U( 50,50) ,V( 50, 50 ) , BETA ( 50 ) , PS I R ( 50) , 

A Ul, VI, DELI ,OELY,RPRES,RDENS,RECOV,ENTP,ND1M,KREF,1REF, 

5 XB(50) ,YB( 50) ,UB(50) ,VB( 50) ,P( 50) ,RECV(2, 10) ,FACTR, I REV, 

6 REG, R( 25,4) , S ( 25 ,4 ) , COWL ( 25 ) , BODY ( 25 ) , NR , NS , ICOWL , I BODY 
DIMENSION AMT( 10) ,THR( 10) , ANG( 10) ,NIS( 10) 

5 READ (5,400) AMO, GAM, TNETC ,BETAE, ALPHA, DELB, DELU1 , ERROR, COWLA, 
1THETAL, START 

READ (5,402) SP I LL , DELP ,OFF , XCOWL , YCOWL , PR I NT 
READ (5,404) NOIM,M, NR ,NS, MTHR ,ND, I SHK, NSHK 
READ (5,406) ( ( R ( I , J ) , J= 1 , 4 ) , COWL ( I ) , I =1 , NR ) 

READ (5,406) < ( S ( I , J ) , J= 1 , 4 ) , BODY ( I ) , I =1 , NS ) 

400 FORMAT ( 6E 12 .0/ 5E 1 2 .0 ) 

402 FORMAT ( 6E 1 2 . 0 ) 
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404 FORMAT (816) 


406 FORMAT (5E12.0) 

408 FORMAT (3E12.0,I6) 

800 FORMAT ( / / 3 7X , 42HCHAR ACT ER I S T I C PROGRAM FOR A BI-RAMP INLET//) 

502 FORMAT ( / / 37X ,42HCH AR ACTER I ST I C PROGRAM FOR A BI-CONE INLET//) 

600 FORMAT ( / / 30X , 3HAM0, IOX , 5HTHETC , 1 OX f 2HBS , I OX , 5HAL PHA , 1 OX , 5HC0WL A, 

\ll) 

602 FORMAT ( 2 3 X , 1 P 5E 1 4 . 5 ) 

604 FORMAT ( / /42 X , 23HC OWL LOCATED AT THETA =lPE12.5//> eLAnru iur 

606 FORMAT ( / / 28X 1 62HSH I FT IN CENTER80DY COEFFICIENTS DOE TO SHOCK IMF 

1INGMENT POINT//) 

608 FORMAT ( / /40X , 6HDELX = 1 PE 1 2 . 5 , 2 X , 6HDELY =IPE12.5//) 

610 FORMAT (1H1) 


REG=1 .0 

SHI FTX = XC OWL -COWL ( 1 ) 
SHIFTY=YCOWL-1.000 
IF ( XCOWL .EQ. 0.0) SHIFTX=0.0 
IF ( YCOWL .EO. 0.0) SHIFTY=0.0 
DO 10 1=1, NR 

CALL SHIFT( I, SHIFTX, SHIFTY) 

10 CONTINUE 

DO 12 K= 1 , 50 


DO 12 1=1,50 

CALL EOUK (K, I ,0.0, 0.0, 0.0, 0.0) 

CALL EOUI ( I ,0.0, 0.0, 0.0, 0.0) 

12 CONTINUE 

CONVP= ( 1 .0+( GAM-1. 0)/2.0)**( -GAM/ (GAM-1.0) ) 
CONVD= ( 1 .0+ (GAM-1 .0) /2.0) #*(-1.0/ ( GAM-1.0) ) 
S I GSQ=( GAM-1.0 )/(GAM+1.0) 

IF (NDIM .GT. 2) GO TO 14 
CALL PSA(AMO,0.0,THETC) 

UC = UK 

VC = VK 

BS=BETAK 

DELC = TAN( THETC ) 

DEL S=TAN ( BS ) 

GO TO 16 

14 CALL CSA(BETAE) 


U1=UC 


V1 = VC 

DEL 1=-1 .O/DELC 
16 CALL OSWR(AMO,BS) 

I R E V= 1 

RECV( 1,IREV)=REC0V 
RECV (2 , IREV ) =RECOV 
F AC T R = R EC 0 V *C 0 N V P 


K = 1 
1 = 1 

S I GN=— 1 .0 
WRITE (6,610) 

IF (NDIM .GT. 2) GO TO 18 
WRITE (6,500) 

GO TO 20 
18 WRITE (6,502) 

20 WRITE (6,600) 

WRITE (6,602) AMO, THETC ,BS , ALPHA, COWL A 
IF (OFF .EO. 0.0) GO TO 26 
XCOWL=COWL ( 1 ) 

YC0WL=1 .000 

THETAL = AT AN (YCOWL /XCOWL ) 
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22 WRITE (6,604) THETAL 
DEL S = DEL P*DELS 
24 CALL OUT 

WRITE (6,610) 

IF ( ABS ( ALPHA ) .GT. .0.0) GO TO 30 
IF (ABS(5(1,2)-S(2,2) ) .GE. ERROR) GO TO 30 
CALL DRAG(THETAL t 1.0,ND) 

DELS=DELP*TAN(BS) 

GO TO 30 

26 IF (SPILL .GT. 0.0) GO TO 28 
THETAL =BS 

28 CALL DRAG(THETAL, 1.0, ND) 

30 K = 1 
1 = 1 

S I GN=— 1 .0 
EXTEND= I SHK 

BODY ( ISHK)=EXTEND*BODY( ISHK) 

R EG=- 1*0 

CALL CONE (K, I , ST ART , ALPHA , S I GN , PR I NT ) 

I REF= 1 

IF (NR .GT. 1 ) GO TO 32 
XS = X ( KR EF , I REF ) 

YS = Y ( KREF , I REF ) 

DELSP=TAN(COWLA) 

XP=XS+3.0 
R EG= 1.0 

CALL CURVE { XS , YS ,DELSP , XP , 0.0, 0.0, 1.0, NR) 

NR=NR+ 1 
GO TO 34 

32 IF ( ABS ( ALPHA ) .GT. 0.0) GO TO 34 
IF (OFF .GT. 0.0) GO TO 34 
XC OWL =C OWL ( 1 ) 

YCOWL = 1 .000 

DIST=XCQWL-X(KREF, IREF) 

IF (ABS(DIST) .GT. ERROR) GO TO 38 
36 CALL EOUMKREF, IREF, XCOWL, YCOWL, UiKREF, IREF) ,V(KREF, IREF) ) 
GO TO 42 

34 IF (OFF .EO. 0.0) GO TO 42 
38 CALL PUT( XCOWL , YCOWL ) 

BODY ( ISHK)=BODY( ISHK) /EXTEND 
42 DELTA=COWLA-ATAN (V (KREF, IREF) /U(KREF, IREF) ) 

IF (ABS(DELTA) .LT. 0.010) DELTA=0.0 
IF (DELTA .GT. 0.0) DELT A=0 .0 
IF ( NSHK .EO. 0) GO TO 64 
IF (NSHK .GT. 1 ) GO TO 44 
S ET P = OF F 
SET Q=OF F 
GO TO 46 
44 SETP=1 .0 
S ET0=0 • 0 
46 CALL MATRX 
KR=KREF 
I R= I REF 

DO 48 J = 1 , NSHK 
S I GN= ( — 1 #0 ) ** ( J+l ) 

REG=S I GN 

IF (NSHK .EQ. 1) BOOY ( I SHK ) =BOOY ( I SHK ) ^EXTEND 
CALL SHOCK(KR, IR, DELTA, SIGN, SETP,SETQ) 

IF (NSHK .EO. 1) BODY( ISHK)=BODY( ISHK)/EXTEND 
KR=KREF 


I R=1 

S TFR = PS I 
DELTA=-DELA 
48 CONTINUE 

IF ( NSHK .GT. 1 ) GO TO 64 
IF (OFF .EO. O.O) GO TO 60 
IF ( I SHK .EQ. 1) GO TO 56 
XSHK=BODY ( I SHK ) 

YSHK = S( ISHK t l )+S( I SHK ,2 )*XSHK+S( I SHK , 3 ) *XSHK**2+S ( I SHK ,4 ) *XSHK**3 

DELX=X(1,1 ) — XSHK 

DELY=Y(1,1)-YSHK 

IF ( ABS ( DEL X ) .GT. ERROR) GO TO 50 
IF ( ABS ( DEL Y ) .LT. ERROR) GO TO 56 
50 R EG=- I . 0 

DO 52 I = I SHK t NS 

CALL SHIFT(I,DELX,DELY) 

52 CONTINUE 

WRITE ( 6*610 ) 

WRITE (6,606) 

WRITE (6,608) DEL X , DEL Y 
54 CALL OUT 

56 IF ( NTHR .EO. 0) GO TO 64 
DO 58 J = 1 , NTHR 
REG=(-1.0)**( J ) 

CALL SHAPE( 1 ,1, REG, 1.0, 1.0) 

58 CONTINUE 
GO TO 66 

60 IF (NTHR .EO. 0) GO TO 64 
KR=KRE F 
I R= I REF 

READ (5,408) ( AMT ( J ) , THR ( J ) , ANG ( J ) , N I S ( J ) , J=1,NTHR) 

DO 62 J= 1 ,NTHR 

CALL SURF(KR,IR,NIS( J) ,AMT( J) , THR ( J ) , ANG ( J ) , STFR , 1 .0 ) 

KR=KREF 
I R = 1 

62 CONTINUE 
64 WRITE (6,610) 

CALL OUT 
66 GO TO 5 
END 
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OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO 


$ I BFTC INLET1 L I ST , REF »DECK » DEBUG 

CHARACTERISTIC PROGRAM FOR ISENTROPIC RAMP INLET, NDIM=2 

CHARACTERISTIC PROGRAM FOR ISENTROPIC SPIKE INLET, NDIM=3 

AMO=FREE STREAM MACH NUMBER 
GAM=RAT 10 OF SPECIFIC HEATS 
THETC=INITIAL CONE OR RAMP ANGLE (RADIANS) 

AMF=EXTERNAL FOCAL POINT MACH NUMBER 
COWLA=COWL ANGLE (RADIANS) 

THETAP=LCOATION OF FOCAL POINT (RADIANS) 

BETAE = EST I MATE OF CONE SHOCK ANGLE (RADIANS) 

R( I ,J)=COEFFICIENTS OF COWL CONTOUR SECTIONS 
S( I ,J)=COEFFICIENTS OF CENTERBODY CONTOUR SECTIONS 
COWL ( I ) =COWL CONTOUR REGIONS 
BODY ( I ) =CENTERBODY CONTOUR REGIONS 

DELB= INITIAL INCREMENT FOR CONE SHOCK ANGLE CALCULATION (RADIANS) 
DELU1=INTEGRATI0N INCREMENT IN CONE FIELD CALCULATION 
ERROR=CONVERGENCE PARAMETER, SET ERR0R=0.1-3 FOR MOST APPLICATIONS 
SH I FTX=SHI FT IN COWL X-COORDINATE FROM DESIGN POSITION 
SH I FT Y = SH I FT IN COWL Y-COORDINATE FROM DESIGN POSITION 

M=NUMBER OF INITIAL I -NET SPACE 
NR=NUM8ER OF COWL CONTOUR SECTIONS 
NS=NUMBER OF CENTERBODY CONTOUR SECTIONS 
NTHR=NUMBER OF ISENTROPIC COMPRESSION SECTIONS 
NSHK=NUM8ER OF INTERNAL SHOCKS 

ND=NUM8ER OF INTEGRATION INCREMENTS IN ADDITIVE DRAG CALCULATION 
I SHK=SHOCK IMPINGMENT POINT ON BODY CONTOUR SECTION 
INDEX=K-INDEX TO INCREASE NUMBER OF GRID MESH POINTS 
NGR I D=NUMBER OF ADDITIONAL GRID SPACINGS 

COMMAND PARAMETERS 

SET NO I M=2 FOR TWU DIMENSIONAL FLOW CALCULATIONS 
SET ND I M=3 FOR THREE DIMENSIONAL FLOW CALCULATIONS 
SET F0CUS=0 . 0 FOR EXTERNAL COMPRESSION TO FOCUS ON SHOCK 
SET FUCUS=1.0 FOR EXTERNAL COMPRESSION TO FOCUS AT THETAP 
SET NSHK=1 FOR ISENTROPIC INTERNAL COMPRESSION 

SET NSHK GREATER THAN 1, ISHK=1, FOR INTERNAL SHOCKS TO REFLECT 
SET 0FF=0 . 0 FOR ON-DESIGN CALCULATIONS 
SET 0FF=1 .0 FOR OFF-DESIGN CALCULATIONS 

COMMON AMO, GAM, THETC, DELB, DELU1 ,DELU,UC, VC, DELC,BS, DELS, ERROR, 

1 0 , AM, AMU, THETA, A,B,C,D,E»F,G,PRES,DENS,AREA,DYNP,CP,PSI, 

2 XK,YK,UK,VK,XJ,YJ,UJ,VJ , X A , YA ,UA , VA , DELA ,BETAK , S IGK , S I GSO , 

3 X( 50,50) ,Y (50,50), U( 50,50) , V f 50, 50 ) , BETA ( 50 ) , PS IR < 50), 

A U1»V1,DEL1,DELY,RPRES»RDENS»REC0V,ENTP,NDIM,KREF,IREF, 

5 XB( 50) ,YB( 50) ,UB( 50) ,VB(50) ,P( 50) ,RECV(2,10) ,FACTR, IREV, 

6 REG, R( 25, A) , S ( 25 , A ) ,COWL ( 25 ) , BODY ( 25 ) , NR ,NS , ICOWL , I BODY 
DIMENSION AMT( 10) ,THR( 10) ,ANG( 10) ,NIS( 10) 

5 READ (5,400) AMO, GAM, THETC ,BE TAE ,AMF , DELB , DELU1 , ERROR ,COWL A , 
1THFTAP, START 

RtAD (5,402) FOCUS, DELP, OFF, XCOWL,YCOWL, PRINT 
READ (5,404) ND IM, M, NR ,NS , NTHR , ND, I SHK , NSHK 
READ (5,406) ( ( R ( I , J ) , J=1 ,4 ) ,COWL ( I ) , I =1 , NR ) 

READ (5,406) ( (S( I ,J ) ,J=1,4) ,BODY ( I ) , 1=1, NS) 

400 FORMAT < 6E12 .0/5E12 .0 ) 

402 FORMAT (6E12.0) 
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404 FORMAT (816) 

406 FORMAT (5E12.0) 

408 FORMAT (3E12.0,I6) 

800 FORMAT ( / / 33X , 50HCMAR ACT ER I ST I C PROGRAM FOR A ISENTROPIC RAMP INLE 
IT//) 

502 FORMAT ( / / 33X , 5 1 HCHAR ACTER I ST IC PROGRAM FOR A ISENTROPIC SPIKE INL 
1 ET / / ) 

600 FORMAT ( / / 30X , 3H AMO, 1 OX , 5 HTHE TC , 1 OX , 2HBS , 1 OX , 5HC0WL A , 1 OX , 3HAMF// ) 
602 FORMAT ( 23X , 1 P5E 14 . 5 ) 

604 F0RMAT(//33X,41HISENTR0PIC COMPRESSION FOCUSES AT THETA =1PE12.5// 
1 ) 

606 FORMAT ( / /28X t 62HSHI FT IN CENTERBODY COEFFICIENTS OUE TO SHOCK IMP 
1INGMENT POINT//) 

608 FORMAT ( //40X , 6HDELX = 1 PE 1 2 . 5 , 2 X , 6HDEL Y =1PE12.5//) 

610 FORMAT ( 1 HI ) 

SHIFTX=XC OWL -COWL ( 1 ) 

SHIFTY=YCOWL— 1 .000 

IF ( XCOWL .EO. 0.0) SH I F TX = 0 • 0 

IF ( YCOWL .EO. 0.0) SHIFTY=0.0 

R EG= 1.0 

DO 10 1=1, NR 

CALL SHIFT ( I, SHIFTX, SHIFTY) 

10 CONTINUE 

DO 12 K = 1 , 50 
DO 12 1=1,50 

CALL EOUK (K, I ,0.0, 0.0, 0.0, 0.0) 

CALL EOUI ( I ,0.0, 0.0, 0.0, 0.0) 

12 CONTINUE 

CONVP= ( 1 .0+ ( GAM-1. 0)/2.0)**( -GAM/ (GAM-1.0) ) 

CONVD=( 1.0+ (GAM-1. 0) /2.0)** (-1.0/ (GAM- 1.0) ) 

S I GSQ= ( GAM-1 .0)/ (GAM+1.0) 

IF (NDIM .GT. 2) GO TO 14 
CALL PS A ( AMO ,0.0,THETC) 

UC = UK 

VC = VK 

BS=BETAK 

DELC = T AN ( THETC ) 

DEL S = T AN ( BS ) 

GO TO 16 

14 CALL CSA ( BETA E ) 

U1=UC 
V1 = VC 

DEL1=— 1 .O/DELC 
16 CALL OS WR ( AMO, BS ) 

I REV= 1 

REC V ( 1, IREV)=RECOV 
RECV(2 ,IREV)=RECOV 
FACTR=RECOV*CONVP 
WRITE (6,610) 

IF (NDIM .GT. 2) GO TO 18 
WRITE (6,500) 

GO TO 20 
18 WRITE (6,502) 

20 WRITE (6,600) 

22 WRITE (6,602) AMO , THE TC , B S , COWL A , AMF 
IF (OFF .EO. 0.0) GO TO 26 
XCOWL=COWL ( 1 ) 

YCOWL=R( 1,1 )+R( 1,2)*XC0WL+R( 1 , 3 ) *XC0WL**2+R ( 1 , 4 ) *XC0WL**3 
24 CALL OUT 

DELQ=DELP*DELS 
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CALL DATUM(START,DELO) 

GU TO 34 

26 IF (FOCUS .GT. O.O) GO TO 28 
THETAP=BS 
GO TO 30 

28 DELS=TAN(THETAP) 

30 WRITE ( 6 , 604 ) THETAP 
K = 1 
1=1 

X(K,I)=1.0/DEL$ 

Y ( K , I ) =1 .0 

CALL DRAG(THETAP,Y(K, I ) ,ND) 

CALL SP IKE ( K , I , M , ST ART ? AMF , 1 .0 , PR INT ) 

I REF=1 

IF (NS .GT. 1) GO TO 32 
XS = X( 1,1 ) 

YS=Y (1,1) 

DELSP=V(1,1 )/U(l,l ) 

XP=XS+3 *0 
R EG=— 1 . 0 

CALL CURVE (XS,YS,DELSP,XP,0. 0,0. 0,1.0, NS) 
NS=NS+ 1 

32 IF (NR .GT. 1) GO TO 34 
XS = X ( KREF * I REF ) 

Y$ = Y ( KREF , I REF ) 

DELSP=TAN( COWL A ) 

XP = XS+3 .0 
R EG= 1.0 

CALL CURVE! XS, YS, DEL SP,XP, 0.0, 0.0, 1.0, NR) 
NR=NR+1 
34 KR= 1 
I R = 1 

REG=-1.0 
EXTENO= I SHK 

BODY ( I SHK ) =BODY ( I SHK) ^EXTEND 
36 CALL SHAPE(KR,IR, -1.0, 0.0, PRINT) 

I R EF = 1 

IF (OFF .EQ. 0.0) GO TO 42 
38 CALL PUT ( XCOWL , YCOWL ) 

42 BODY ( I S HK ) =BOD Y ( ISHK)/EXTEND 

OELTA=COWLA-ATAN (V (KREF, I REF) /U ( KREF, IREF) ) 
IF ( ABS (DELTA) .LT. 0.010) DEL T A = 0 . 0 
IF (DELTA .GT. 0.0) DELTA=0 .0 
IF ( NSHK .EO. 0) GO TO 64 
IF (NSHK .GT. 1) GO TO 44 
S ETP=OFF 
SET 0=0FF 
GO TO 46 
44 $ETP=1 .0 
S ET0 = 0 • 0 
46 CALL MATRIX 
KR=KREF 
I R= I R EF 

DO 48 J«1,NSHK 
SIGN=(-l .0)**( J+l ) 

REG=S I GN 

CALL SHOCK (KR,IR, DELTA, SI GN,SETP,SETQ) 

KR=KREF 

IR = 1 

S TFR = PS I 


DELTA=-DELA 
48 CONTINUE 

IF ( NSHK .GT. 1) GO TO 64 
IF (OFF .EG. O.O) GO TO 60 
IF ( ISHK .EO. 1 ) GO TO 56 
XSHK=BODY ( ISHK) 

YSHK=S ( I SHK » 1 )+S( I SHK ,2 )*XSHK+S( I SHK , 3 )»XSHK**2+S ( I SHK ,4 ) *XSHK**3 

DELX=X( 1,1 )-XSHK 

DELY=Y(1,1)-YSHK 

IF (ABS(DELX) .GT. ERROR) GO TO 50 
IF ( ABS ( DEL Y ) .LT. ERROR) GO TO 56 
50 REG=— 1.0 

DO 52 I=ISHK,NS 

CALL SHIFT( I ,DELX,DELY) 

52 CONTINUE 

WRITE (6,610) 

WRITE (6,606) 

WRITE (6,608) DELX, DEL Y 
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CALL OUT 



IF 

(DELX , 

.LT. 0.0) GO TO 56 


GO 

TO 66 


56 

IF 

(NTHR , 

.EO. 0) GO TO 64 


DO 

58 J = 1 

,NTHR 


REG 

= <-1.0 

) ** ( J ) 


CALL SHAPE(1,1,REG,1.0,1.0) 

58 CONTINUE 
GO TO 66 

60 IF ( NTHR .EO. 0) GO TO 64 
KR=KREF 
I R = I REF 

READ (5,408) ( AMT ( J ) , THR ( J ) , ANG ( J ) ,NI S ( J ) , J=1,NTHR) 

DO 62 J=1 ,NTHR 

CALL SURF(KR»IR»NIS(J) , AMT ( J ) , THR ( J ) , ANG ( J ) , STFR . 1 .0 ) 
KR=KREF 
I R = 1 

62 CONTINUE 
64 WRITE (6,610) 

CALL OUT 
66 GO TO 5 
END 
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SIBFTC EQUKS L I ST , REF , DECK , DEBUG 

C 

C EQUIVALENCE SUBROUTINE 

C 

SUBROUTINE EQUK ( KR , I R , XA , YA , UA , VA ) 

COMMON AMO,GAM,THETC,DELB,DELUl,DELU,UC,VC,DELC,BS, DELS, ERROR, 

1 Q,AM,AMU,THETA,A,B,C,D,E,F,G,PRES,DENS,AREA,DYNP,CP,PSI , 

2 XK,YK,UK,VK,XJ, YJ,UJ, VJ,XA, YA,UA, VA, DELA , BETAK , S IGK , S I GSQ, 

3 X ( 50 * 50 ) ,Y (50,50 ) ,U( 50,50) f V( 50,50) ,BETA< 50) ,PSIR( 50) , 

4 U1 ,VI,DELI,DELY,RPRES,RDENS,RECOV,ENTP t NDIM f KREF, I REF, 

5 XB(50),YB(50),UB(50) ,VB(50) ,P(50) ,RECV(2,10) r FACTR r IREV, 

6 REG, R( 25,4) , S ( 2 5 ,4 ) , COWL ( 2 5 ) ,B0DY<25) ,NR,NS, I COWL t I BODY 
X ( KR , I R ) =X A 

Y ( KR , I R ) =Y A 
U ( KR , I R ) =UA 
V ( KR , I R ) =VA 

IF ( X ( KR , I R > .EQ. 0.0) GO TO 10 

CALL CFPR(X(KR,IR),Y(KR,IR),U(KR, IR),V(KR,IR) ) 

10 RETURN 
END 


$ I BFTC EOUIS LIST, REF, DECK, DEBUG 

C 

C EQUIVALENCE SUBROUTINE 

C 

SUBROUTINE EQU I ( I R , X A , Y A ,UA, V A ) 

COMMON AMO, GAM, THETC, DELB, DELU1, DEL U,UC, VC, DELC,BS, DELS, ERROR, 

1 Q,AM,AMU,THETA,A,B,C,D,E,F,G,PRES,DENS,ARFA,DYNP,CP,PSI , 

2 XK,YK,UK,VK,XJ»YJ,UJ»VJ,XA,YA,UA,VA,DELA,BETAK,SIGK,SIGSQ, 

3 X(50,50) ,Y(50,50) ,U(50,50) ,V(50,50) ,BETA(50) ,PSIR(50) , 

4 U1 , VI, DELI ,DELY,RPRES,RDENS,RECOV,ENTP,NDIM,KREE, I REF, 

5 XB ( 50) ,YB(50) ,UB(50) ,VB( 50) ,P(50) ,RECV(2, 10) , FACTR , I REV , 

6 REG»R(25,4) ,S(25,4) ,C0WL(25) ,BODY (25) ,NR,NS, ICOWL, I BODY 
XB ( I R ) =XA 

YB ( I R ) =YA 
UB ( I R ) =UA 
VB ( I R ) =V A 

IF ( XB ( IR) .EQ. 0.0) GO TO 10 

CALL CFPR( XB< IR ) ,YB( IR ) ,UB< IR) ,VB( IR) ) 

10 RETURN 
END 
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$ I BFTC CURVES L I ST , REF , DECK , DEBUG 
C 

C ROUTINE FOR CALCULATING COEFFICIENTS OF SURFACE CONTOURS 

C 

SUBROUTINE CURVE ( X 1 , Y 1 , SL PI , X2 , Y2 , S L P2«, EX P , I NDEX ) 

COMMON AMO, GAM, THETC, DEL B, DE LU1, DELU, UC, VC, DELC,BS, DELS, ERROR, 

1 0 , AM , AMU , THETA, A, B,C,D,E,F,G, PRES, DENS, AREA, DYNP, CP, PSI, 

2 XK,YK,UK,VK,XJ,YJ,UJ,VJ,XA,YA,UA,VA,DELA,BETAK,SIGK,SIGSO, 

3 X( 50,50) ,Y<50,50) ,U( 50,50) ,V( 50,50 ) ,BETA( 50) ,PSIR( 50) , 

4 U I, VI, DEL I ,DELY,RPRES,RDENS,RECOV,ENTP,NDIM,KREF, I REF, 

5 XB ( 50 ) , YB ( 50 ) , UB ( 50 ) , V B ( 50 ) , P ( 50 ) , R EC V ( 2 , 1 0 ) , F ACTR , I RE V , 

6 REG, R( 25,4) ,S(25,4) ,COWL(25) ,B0DY(25) ,NR,NS, I COWL, I BODY 
Z = X2 —X I 

IF (Z .GT. 0.0) GO TO 8 
A0 = 0 . 0 
A 1 =0 • 0 
A2=0 .0 
A 3 = 0.0 
GO TO 14 
8 A0 = Y1 
A1 = SL P 1 
DYR=SLP2 

IF (EXP .EO. 3.0) GO TO 10 
IF (EXP .EO. 2.0) GO TO 12 
A2 = 0 . 0 
A 3 = 0 • 0 
GO TO 14 

10 FUN1=Y2-Y1-A1*Z 
FUN2=DYR-A1 
DEL T A=Z*#4 

A2=( 3,0*FUN1-Z*FUN2 )*Z**2/ DELTA 
A3=( Z*FUN2-2.0*FUN1 )*Z /DELTA 
GO TO 14 

12 A2=(DYR-A1 )/(2.0*Z> 

A3=0 .0 

14 P0=A0-A1*X1+A2*X1**2-A3*X1**3 
P1=A1-2.0*A2*X1+3.0*A3*X1**2 
P2=A2-3.0*A3*X1 
P 3 = A 3 

IF (REG .EO. “1.0) GO TO 18 
R ( INDEX, 1)=P0 
R( INDEX, 2 ) =P 1 
R ( INDEX, 3)=P2 
R ( INDEX, 4)=P3 
COWL ( I NDEX ) =X 1 
C OWL ( INDEX+1 ) =X2 
DO 16 J = 1 , 4 
R ( INDEX+1 , J >=0.0 
16 CONTINUE 
GO TO 22 

18 S( INDEX, 1)=P0 
S( INDEX , 2 ) =P 1 
S( INDEX, 3)=P2 
S< INDEX, 4)=P3 
BODY ( I NDEX ) = X 1 
BODY ( INDEX+1 )=X2 
DO 20 J=l,4 
S( INDEX+1, J)=0.0 
20 CONTINUE 
22 RETURN 
END 



$ I B FTC SHIPS LIST, REF, DECK, DEBUG 

C 
C 
C 

SUBROUTINE SHI FT ( I , DEL X,DELY ) 

COMMON AMO,GAM,THETC,DELB,DELUl ,DELU,UC, VC, DELC,BS, DELS, ERROR, 

1 0, AM, AMU, THETA, A, B,C,D,E,F,G, PRES, DENS, AREA, DYNP, CP, PSI, 

2 XK , YK ,UK , VK ,* J , Y J , U J , V J , X A , YA , UA , VA , DEL A , BET AK , S I GK , S I GSO , 

3 X ( 50,50) ,Y( 50,50) ,U( 50,50) , V ( 50 , 50 ) , BETA ( 50 ) , PS I R ( 50 ) , 

A U1 ,V1 , DEL 1, DEL Y, RPRE S,R DENS, RECOV, ENT P,NDIM,KREF, I REF, 

5 XB ( 50) ,YB( 50) , UB ( 50 ) , VB < 50 ) , P ( 50 ) , RECV ( 2 , 1 0 ) , FACTR , I REV, 

6 REG, R (25,4) , S ( 25 ,4 ) , COWL! 25 ) , BODY ( 2 5 ) , NR , NS , I COWL , I BODY 
IF (REG • EQ • -1.0) GO TO 10 

P1=R( I , 1 ) +DEL Y-R ( 1,2) *DELX+R ( I , 3 ) *DELX**2-R ( I ,4)*DELX**3 
P2 = R ( I ,2 )-2.0*R( I ,3)*DELX+3.0*R( I ,4)*DELX**2 
P3 = R( I ,3)-3.0*R( I , A ) *DEL X 
P4=R( 1,4) 

R( I ,1 )=P1 
R ( 1 , 2 ) = P2 
R ( I , 3 ) =P3 
R( I ,4) =P4 

COWH I ) =COWL ( I ) +DEL X 
GO TO 12 

10 P1=S( 1,1) +DEL Y-S ( 1,2) ^DELX+S ( 1,3) #DELX**2-S ( I , 4 ) #DELX**3 
P2=S ( I ,2 )-2.0*S ( 1,3) *DELX+3. 0#S ( 1,4) *DELX**2 
P3=S ( I ,3)-3.0*S( I ,4)*DELX 
P4=S ( 1,4) 

S( I ,1 )=P1 
S( I , 2 ) = P2 
S( I , 3 ) =P3 
S ( I , 4 ) = P4 

BODY ( I ) =B ODY ( I )+DELX 
12 RETURN 
END 


$ I BFTC OUTR LIST, REF, DECK, DEBUG 


C 

C ROUTINE FOR WRITING OUT TABLE OF CONTOUR COEFFICIENTS 

C 


100 

102 

104 

106 

108 


10 


12 


1 

2 

3 

4 

5 

6 


SUBROUTINE OUT 

COMMON AMO, GAM, THETC, DELB, DELU 1 ,DELU,UC, VC, DELC,BS, DELS, ERROR, 
0,AM,AMU,THETA,A,B,C,D,E,F,G,PRES,DENS,AREA,DYNP,CP,PSI, 
XK,YK,UK,VK,XJ,YJ,UJ,VJ,XA,YA,UA,VA,DELA,BETAK,SIGK,SIGSQ 
X( 50,50) ,Y( 50,50) ,U(50, 50) ,V(50,50) ,BETA( 50) ,PSIR( 50) , 

U1 ,V1 ,DEL1,DELY,RPRES,RDENS,REC0V,ENTP,NDIM,KREF, IREF, 
XB(50),YB(50) ,UB(50),VB(50) , P ( 50 ) , R ECV ( 2 , 1 0 ) , F ACTR , I REV , 
REG,R(25,4) , S ( 25 ,4 ) , COWL ( 25 ) ,BODY(25) , NR , NS , I COWL , I BODY 
FORMAT (//45X,26HTABLE OF COWL COEFFICIENTS//) 

FORMAT ( //30X,2HR1,12X,2HR2,12X,2HR3,12X,2HR4,11X,4HC0WL// ) 
FORMAT ( 23X , 1 P5E 14 .5// ) 

FORMAT ( //43X,32HTABLE OF CENTERBODY COEFFICIENTS//) 

FORMAT ( / / 30X , 2HS 1 , 12X,2HS2,12X,2HS3, 12X, 2HS4, 10X, 4HB0DY// ) 

WRITE (6,100) 

WRITE (6,102) 

DO .10 1=1, NR 

WRITE (6,104) (R( I,K) ,K = 1,4) ,COWL( I ) 

CONTINUE 
WRITE (6,106) 

WRITE (6,108) 

DO 12 1=1, NS 

WRITE (6,104) (S(I,K),K=1,4),B0DY( I ) 

CONTINUE 

RETURN 

END 


T 
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$ I BFTC FINDS L I ST , REF , DECK , DEBUG 

C 
C 
C 

SUBROUTINE FIND(POINT) 

COMMON AMO,GAM,THETC,DELB,OELUl,DELU,UC, VC, DELC,BS, DELS, ERROR, 

1 O t AM,AMU,THETA,A,B,C,D,E,F,G,PRES,DENS,AREA,DYNP,CP,PSI , 

2 XK,YK,UK,VK,XJ,YJ,UJ,VJ,XA,YA,UA,VA,DELA,BETAK,SIGK,SIGSQ» 

3 X< 50,50) ,Y( 50,50 >,U( 50, 50 ) , V ( 50, 50 ) , BETA ( 50) ,PSIR( 50), 

4 U1 , VI, DELI ,DELY,RPRES,RDENS,RECOV,ENTP,NDIM,KREF, I REF, 

5 XB(50),YB(50)*UB(50),VB(50)»P(50),RECV(2,10),FACTR,IREV, 

6 REG,R ( 25,4 ) ,S( 25,4) ,COWL( 25 ) , BODY { 25 ) , NR , NS , I COWL, I BODY 
100 FORMAT ( //25X,52HC0WL POINT LIES OUTSIDE RANGE OF INPUT CONTOURS, 

IX =1PE12.5//) 

102 FORMAT ( / /2 5X , 52HB0DY POINT LIES OUTSIDE RANGE OF INPUT CONTOURS, 
IX =1 PE 12 . 5// ) 

MR=NR— 1 
MS=NS-1 

IF (REG .EO. -1.0) GO TO 18 

REF=COWL ( 1 l-POINT 

IF (REF .GT. 0.0) GO TO 14 

DO 10 1=1, MR 

REF=COWL( 1+1 )-POINT 



I F 

(REF .GT. 

0.0) 

GO 

TO 

12 

10 

CONTINUE 





12 

IF 

(REF ,01. 

0.0) 

GO 

TO 

16 

14 

WRITE (6,100) 

POINT 




CALL CFPR (0.0, 0.0, 0.0, 0.0) 
CALL EXIT 
16 ICOWL= I 
GO TO 28 

18 REF=BODY( 1 l-POINT 



IF (REF .GT . 
DO 20 1=1, MS 

0.0 ) 

GO 

TO 

24 


REF=BODY ( 1+1) 

-POINT 




IF (REF .GT. 

0.0) 

GO 

TO 

22 

20 

CONTINUE 





22 

IF (REF .Gf. 

0.0) 

GO 

TO 

26 

24 

WRITE (6,102) 

POINT 




CALL CF PR ( 0 . 0 
CALL EXIT 

,0.0 

,0. 

• 

o 

► 

o 

0) 

26 

I BODY= I 





28 

RETURN 






END 
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I 


$ I BFTC CCRAS L I ST , REF , DECK , DEBUG 

C 

C SUBROUTINE FOR STARTING SHOCK SOLUTION 

C 

SUBROUTINE CC R A ( KQ t I Q t XP t Y P t UP t V P , A NGL E f S I GN > 

COMMON AMOtGAM, THETC» DELB tDELUItDELU, UC t VC t DELC»BSt DELS »ERROR t 

1 OtAMfAMUt THETA t A tB tC » Dt E t F tG t PRES t DENS t AREA tDYNP tCPtPSI t 

2 XK f YK f UK f VK t XJ f YJtUJ-t V J , X A , YA , UA , V A , DEL A t BE T AK t S I GK t S I GSQ* 

3 X< 50,50) ,Y< 50,50) ,U( 50,50) * V ( 50 , 50 ) , BETA ( 50 ) t PSIR(50) , 

4 U1 , VI ,UEL1 ,DELY,RPRES,RDENS,RECOV,ENTP,NDIM,KREF, I REF , 

5 XB< 50) ,YB(50) ,UB(50) ,VB(50) ,P(50) ,RECV(2,10) ,FACTR,IREV, 

6 REG,R<25,4) ,S( 25,4) , COWL (25) , BODY ( 25 ) f NR , NS , I COWL, I BODY 

100 FORMAT (//39X,47H UNABLE TO OBTAIN CONVERGENCE IN SUBROUTINE CCRA// 
1 ) 

102 FORMAT ( / / 32 X f 4HXR EF , 10X , 4H YR EF , 10X , 3HAMR , 1 1 X f 4HQR EF t 9X , 6HTHET AR / / 
1 ) 

104 FORMAT ( 2 8 X , 1 P 5E 1 4 . 5 / / ) 

XR E F = X P 
YR E F = YP 
DELTA=ANGLE 

CALL CFPR { XP,YP,UP,VP) 

AMP= AM 

THET A P = THET A 
OR E F = 0 
I T ER = 0 

10 T EST = OR EF 
I T ER= I T ER+ 1 

CALL PSA< AMP, THETAP, DELTA) 

UR EF = UK 
VR EF = VK 

12 CALL CCRE(XREF,YREF ,UREF , VREF , S I GN ) 

X J = XK 
YJ = YK 
U J = UK 
V J = VK 

IF (ABS(ANGLE) . LT . 0.001) GO TO 22 
14 CALL CCRE(XREF,YREF ,UREF , VREF , -S IGN ) 

XQ = XK 
YQ = YK 

TANQ=VK/UK 
DU=(UJ-U(KG, 10) ) 

DX=(XJ-X(K0,IQ) ) 

DY=( YJ-Y(KQtlO) ) 

DS=SORT ( DX**2+DY**2 ) 

DUDS=DU/DS 

DELP=XQ-X(KQ,IQ) 

DELO= YO-Y ( KQ , I Q ) 

DELR = SQRT ( DEL P**2+DELQ**2 ) 

UQ=U(KO,IG)+DUDS*DELR 

VQ=TANQ*UQ 

16 CALL CCRC(XJ,YJ,UJ,VJ,XQ,YG,UQ,VQ,SIGN) 

XRE F = X P 
YREF=YP 
URFP=UK 
VREF=VK 

18 CALL CFPR(XREF,YREF,UREF,VREF) 

0REF=0 

DELTA=THETA-THETAP 

IF ( ABS ( DELTA ) .LT. ERROR) DELT A=0 .0 
IF ( ITER .LT. 25) GO TO 20 
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WRITE (6,100) 

WRITE (6,102) 

CALL CFPR( XREF,YREE,UREF,VREF) 

WRITE (6,104) XREF,YREF, AM, Q, THETA 
CALL EXIT 

20 IF ( ABS ( TEST-ORE F ) ,.GT. ERROR*QREF) GO TO 10 
22 XK=XREF 
YK=YREF 
UK=UREF 
VK= VREF 
RETURN 
END 


SIBFTC CCRBS 
C 

C SUBROUT 


LIST, REF, DECK, DEBUG 


SUBROUT 

COMMON 

1 

2 

3 

4 

5 

6 

100 FORMAT 
1) 

102 FORMAT 
1 ) 

104 FORMAT 
XR E F= X P 
YREF=YP 



( //34X,4HXREF,10X,4HYREF,10X,3HAMR,11X,4H0REF,9X,6HTHETAR// 
( 28X , 1 P5E 14*5//) 


IP=IQ 

DELTA = ANGL E 

CALL CFPR(XP,YP,UPrVP) 

AMP=AM 

THET AP = THET A 
OR EF = 0 
I TER=0 

10 TEST=OREF 
I TER= I TER+1 

CALL PSA(AMP,THETAP, DELTA) 

UREF-UK 
VREF= VK 

CALL CFPR(XREF,S I GN*YRE F , UR EF , S I GN#VREF ) 

SLOPE=SIGN*B 

12 CALL CCRC(XREF,YREF,UREF,VREF,X(KO, 10) ,Y (KO, 10) ,U(KO, 10) , 
1 V(KOtlO) ,SIGN) 

X J = XK 

Y J = YK 
U J - UK 

V J = VK 

IF ( ABS ( ANGLE ) .LT. 0.001) GO TO 22 
IF ((XJ-XREF) .LT. 0.0) GO TO 22 
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14 CALL BUUND(KQ,1 ,KQ, I P , XREF , YREF , SLOPE, SIGN) 
XG=XK 
YO = YK 
UO=UK 
V0 = VK 

REF=X (KO, I P )-X0 

IF (REF . GE • O.O) GO TO 16 

IP=IP+1 

GO TO 14 

16 CALL CCRC(XJ,YJ,UJ,VJ,XO,YO,UQ,VQ,SIGN) 

XR E F = X P 
YR EF = Y P 
UREF=UK 
VR EF = VK 

18 CALL CFPR(XREF, YREF , UREF , VREF ) 

OR EF = Q 

DELTA=THETA-THETAP 

IF (ABS(DELTA) .LT. ERROR) DELT A=0 .0 
IF ( ITER .LT. 100) GO TO 20 
WRITE (6,100) 

WRITE (6,102) 

CALL CFPR( XREF, YREF,UREF, VREF) 

WRITE (6,104) XREF, YREF, AM, 0, THETA 
CALL OUT 

CALL CFPR(0. 0,0. 0,0. 0,0.0) 

CALL EXIT 

20 IF ( ABS ( TEST-QREF ) .GT. ERROR*QREF) GO TO 10 
22 XK = XR E F 
YK = YR E F 
UK=UREF 
VK = VR E F 
RETURN 
END 


$ I BFTC CCRCS L 1ST, REF, DECK, DEBUG 


C 

c 

c 


CONIC CHARACTERISTIC ROUTINE C 



Y I =S I GN# Y2 
UI=U2 

V I =S I GN*V2 
X F = X3 

YF=SIGN*Y3 
UF = U3 

VF=SIGN*V3 
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8 CALL C F PR (XI,YI,UI,VI ) 

A I = A 

B I = B 

D I = D 

El = E 

FI = F 

G I = G 

01=0 

CALL CFPR(XF,YF,UF,VF) 

AJ = A 

BJ = B 

0J = D 

EJ = E 

FJ = F 

GJ = G 

0J = 0 

AK = .5*(AI + AJ) 

BK = .5*(BI + 8 J ) 

DK = .5*(DI + OJ) 

EK = .5*(EI + EJ) 

FK = .5*(FI + FJ) 

GK = .5*<GI + GJ) 

0K= • 5* ( 0 I+Q J ) 

10 I T ER= 1+ I T ER 

A = .5 * ( A I + AK) 

B = .5 * (BJ + BK) 

D = • 5 * ( D I + DK) 

E = . 5 * ( E J + EK) 

F = •5* (FI + FK) 

G = .5 * (GJ + GK) 

XK= ( (YF-YI )+A*XI-B*XF)/(A-B) 

YK= YF+B*(XK-XF) 

VK= ( (UF-UI )+D*VI-E*VF+G*(XK-XF)-F* ( XK-X I ) ) / ( D-E ) 
UK= UF+E*(VK-VF)+G*(XK-XF) 

TEST = 0K 

CALL CFPR (XK,YK,UK,VK) 

AK = A 
BK = B 
DK = D 
EK = E 
FK = F 
GK = G 
0K = 0 

IF ( ITER .LT. 25) GO TO 12 
WRITE (6,100) 

CALL CFPR(0. 0,0*0, 0.0, 0.0) 

12 IF ( ABS ( TEST-OK ) .GT. ERROR*OK) GO TO 10 
XK = XK 

YK=SIGN*YK 
UK = UK 

VK=SIGN*VK 
14 RETURN 
END 
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LIST, REF, DEBUG 


$ I BFTC CORES 
C 

C CHARACTERISTIC ROUTINE FOR A CURVED SURFACE 

C 

SUBROUTINE CC RE ( X3 , Y3 , U3 , V3 , S I GN ) 

COMMON AMO, GAM, THETC, DELB, DELUI, DELU,UC, VC, DELC,BS, DELS, ERROR, 

1 0,AM,AMU,THETA,A,B,C,D,E,F,G,PRES,DENS,AREA,DYNP,CP,PSI , 

2 XK,YK,UK,VK,XJ,YJ,UJ, V J , XA , YA, UA, VA, DELA , BETAK, S IGK , S IGSQ, 

3 X( 50,50) ,Y (50,50 ) ,U( 50,50) ,V( 50,50) ,BETA( 50) ,PSIR( 50) , 

A Ul, VI, DELI ,DELY,RPRES,RDENS,RECOV,ENTP,NDIM,KREF, I REF, 

5 XB(50),YB(50),UB(50),VB(50) , P ( 50 ) , REC V ( 2 , 10 ) ,F ACTR , I REV , 

6 REG, R( 25,4) ,S<25,4) ,C0WL(25) ,BODY (25) ,NR,NS, ICOWL, I BODY 
100 FORMAT ( //39X,38HIMPR0PER CALCULATION PERFORMED IN CCRE//) 

102 FORMAT ( //22X , 2HP 1 , 12X , 2HP2 , 12X, 2HP3, 12 X , 2HP4, 1 2X , 2HP5 , 12 X, 2HP6// ) 
104 FORMAT ( 1 5X , 1 P6E 14 .5// ) 

106 FORMAT ( //48X, 7HP0INT =1PE12.5//) 

108 FORMAT ( / / 39 X , 47HUN ABL E TO OBTAIN CONVERGENCE IN SUBROUTINE CCRE// 
1) 

I TER = 0 
XREF=X3 
XC = X3 

YC=S I GN*Y3 
UC=U3 

VC=S I GN*V3 

CALL CFPR(XC,YC,UC,VC) 

AK = A 
DK = D 
FK = F 
AR EF = AK 
UREF-UC 
VREF=VC 
0REF=0 
POI NT = X3 
10 I TER= 1+ I TER 

CALL FIND(XREF) 

IF (REG .GT. 0.0) GO TO 12 
AO = SIGN*S( I BODY , 1 ) 

A1 = SIGN*S( I BODY , 2 ) 

A2 = SIGN*S( I BODY , 3 ) 

A3=SIGN*S< I BODY ,4 ) 

GO TO 14 

12 A0= SIGNER ( ICOWL ,1 ) 

A1=SIGN*R( ICOWL, 2) 

A2 = S I GN*R ( ICOWL, 3) 

A3=S I GN*R ( ICOWL, 4) 

14 B0=A0— YC+AREF&XC 
B1=A1-AREF 

CALL R00T3(A3,A2,B1,B0) 

M=0 

DO 16 1=1,3 
J=2*I-1 
K = 2* I 

P ( J ) = P ( J ) -PO I NT 

IF ( ABS ( P ( K ) ) .GT. ERROR) GO TO 16 
M=M+ 1 

PSIR(M)=P(J) 

16 CONTINUE 
XMIN=10.0 

IF (M .GT. 0) GO TO 18 
WRITE (6,100) 
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18 


WRITE (6,106) POINT 
WRITE (6,102) 

WRITE (6,104) ( P ( I ) , I =1 ,6 ) 

CALL CFPR(0. 0,0. 0,0. 0,0.0) 

CALL EXIT 
DO 20 1=1, M 
DI ST=ABS ( PS I R ( I ) ) 

XM I N=AM INI ( XM IN, 01 ST ) 

20 CONTINUE 
DO 22 1=1, M 
DI ST =ABS ( PS IP ( I ) ) 

IF (DIST .EO. XMIN) GO TO 24 
22 CONTINUE 

24 XREF=PS I R ( I (+POINT 
YREF=A0+A1*XREF+A2*XREF*#2+A3*XREF#*3 
CALL CFPR( XREF,YREF,UREF,VREF ) 

TEST=0 

DREF=0.5*(DK+U) 

FREF=0.5*(FK+F) 

SL0PE=A1+2.0*A2*XREF+3.0*A3*XREF**2 

UREF= (UC-DREF*VC+FREF*( XREF-XC ) ) / ( 1 .0-DREF*SL0PE ) 

VREF=UREF*SLOPE 

CALL CFPR( XREF,YREF,UREF , VREF ) 

0REF=0 

AREF=0.5*(AK+A) 

IF ( ITER .LT. 25) GO TO 25 
WRITE (6,108) 

CALL CFPR(0. 0,0. 0,0. 0,0.0) 

25 IF ( ABS ( TEST-ORE F ) .GT. ERROR*GREF) GO TO 10 
XK=XREF 

YK=SIGN#YREF 

UK=UREF 

VK = S I GN*VRE F 

26 RETURN 
END 


$ I BFTC MATRIX L I S T , RE F , DEC K , DE BUG 
MATRIX INVERSION SUBROUTINE 
SUBROUTINE MATRX 

COMMON AMO, GAM, THETC, DEL B, DELU1 ,DELU,UC, VC, DELC,BS, DELS, ERROR, 

1 0, AM, AMU, THETA, A, B , C , D, E , F , G , PR ES , DENS , ARE A, DYNP , CP , PS I , 

2 XK,YK,UK,VK,XJ,YJ,UJ,VJ»XA,YA,UA,VA,DELA,BETAK,SIGK,SIGSQ, 

3 X( 50, 50) ,Y( 50,50) ,U( 50,50) ,V( 50,50) ,BETA( 50) ,PSIR( 50) , 

4 U1,V1 ,DEL1,DELY,RPRES,RDENS,REC0V,ENTP,NDIM,KREF, I REF, 

5 XB(50),YB(50),UB(50),VB(50),P(50) , R ECV ( 2 , 1 0 ) , FACTR , I REV , 

6 REG»R (25,4) , S ( 25 ,4 ) , COWL ( 25 ) , BODY (25) , NR, NS , I COWL , I BODY 
DO 12 K= 1 , 50 
DO 10 1=1,50 
IF (I .LE. K) GO TO 10 
CALL EOUK (K,I ,0.0, 0.0, 0.0, 0.0) 

10 CONTINUE 
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12 CONTINUE 

DU 16 K=l*50 
UO 14 1=1*50 
IF (I .EQ. K) GO TO 16 

CALL EOUK ( I *K*X<K*I ) *Y(K* I ) *U(K, I ) ,V(K* I ) ) 
14 CONTINUE 
16 CONTINUE 

00 20 1=1*50 

DO 18 K = 1 * 50 

IF (K .LE. I ) GO TO 18 

CALL EGUK(K*I *0.0*0. 0*0. 0*0.0) 

18 CONTINUE 
20 CONTINUE 
RETURN 
END 


$ I BFTC CORP LIST*REF*DECK* DEBUG 
C 

C CALCULATION OF LOCAL STREAMLINE CONDITIONS 

C 

SUBROUTINE CORE ( K * I * S I GN ) 

COMMON AMO, GAM *THETC* DELB *DELU1 *DELU*UC * VC, DELC*BS* DELS* ERROR* 

1 0*AM*AMU*THETA,A*B,C*D*E*F*G*PRES*DENS*AREA*DYNP*CP*PSI* 

2 XK * YK *UK * VK * X J * Y J * U J * V J * X A * Y A * UA * VA * DEL A * BET AK * S I GK * S I GSO 

3 X(50*50),Y(50*50) ,U( 50*50) * V ( 50 * 50 ) * BETA ( 50 ) * PS I R ( 50 ) * 

4 U1 *V1 *DEL1 * DEL Y, RPRES, R DENS *RECOV* ENT P*NDIM*KREF* I REF* 

5 XB(50) *YB(50) *UB(50) *VB(50) * P ( 50 ) *RECV(2*10) * F ACT R * I REV * 

6 REG *R( 25*4) * S ( 25 *4 ) * COWL ( 25 ) * BODY ( 2 5 ) . NR * NS * I COWL * I BODY 
CALL CFPR(X(K* I )*Y(K* I ) *U(K* I ) ,V(K*I) ) 

. PRHO=PRES*FACTR 
THET AR = THET A 
10 FACTA=(GAM~1.0)/2.0 
F ACTB= ( GAM+1 .0 ) /2 • 0 
EXP=- ( G AM— 1 .0 ) /GAM 
CON V P= ( 1 .0+FACTA)**< 1.0/ EXP) 

RAVE = FACT R /C ONVP 
K S I GN= S I GN 
KREV=(3+KSIGN)/2 
RATIO=RAVE/RECV(KREV* IREV) 

IF (ABS(RATIO-l.OOO) .LE. ERROR) GO TO 14 
PRH=( PRHO/RAVE)*RATIO 
AMRSO= ( (PRH)**EXP-1.0)/FACTA 
OR SO=F AC TB^AMRSO / ( 1 . 0+F AC TA* AMRSO ) 

OR=SORT ( ORSO ) 

UR=QR*COS(THETAR) 

VR=QR*SIN( THETAR ) 

CALL EOUK (K*I*X(K*I ) * Y ( K * I ) * UR * VR ) 

IF (I .LT. K) GO TO 14 
12 CALL EOUI ( I*XB( I ),YB( I )*UR*VR) 

14 RETURN 
END 
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$ I BFTC ROOTE L I ST t REF t DECK , DEBUG 
C 

C DETERMINATION OF THE ROOTS OF A CUBIC EQUATION 

C 

SUBROUTINE R00T3 ( Z I , T2 , l 3 , Z4 ) 

COMMON AMO,GAM,THETC,DELB,DELUl, DELU, UC , VC , DELC , BS , DELS, ERROR, 

1 Q, AM, AMU, THETA, A , B , C , D , E , E , G , PR E S , D ENS , AR E A , DYNP , CP, PS I , 

2 XK,YK,UK,VK,XJ, YJ,UJ,VJ,XA, Y A , UA , V A , DEL A , BE TAK , S I GK , S I GSQ, 

3 X( 50,50) ,Y (50,50) ,U( 50,50) ,V(50,50) , BETA* 50) , PS IR (50) , 

4 U1 ,V1 ,DEL1 ,DELY,RPRES,RDENS,RECOV,ENTP,NDIM,KREF, I REE, 

5 XB(50) ,YB(50) ,UB( 50) ,VB(50) ,P(50) ,RECV(2,10) ,FACTR,IREV, 

6 REG,R(25,4),S<25,4) , COWL ( 2 5 ) , BODY ( 25 ) , NR , NS , ICOWL, I BODY 
DIMENSION AUG ( 10) t AL PH I ( 10) , BET I ( 10) ,RG( 10) 

100 FORMAT (//43X,39HIMPR0PER IMPUT DATA TO SUBROUTINE R00T3//) 

102 FORMAT ( / /42 X , 2HZ 1 , 1 2 X , 2HZ 2 , 12 X , 2HZ 3 , 1 2X , 2HZ4/ / ) 

104 FORMAT ( 34X * 1 P4E 14. 5// ) 

IF (ABS(Zl) .LT . 0.00001) GO TO 26 
R2 = Z 2/ Z 1 
R3 = Z 3/ Z 1 
R4=Z4/Z1 

A=(3.0*R3-R2**2)/3.0 

B=( 2.0*R2**3-9.0*R2*R3+27.0*R4) /27.0 
TEST=B**2/4.0+A**3/27.0 
EXAM =-B/2.0 
CHECK=SQRT( ABS( TEST ) ) 

IF (CHECK .GT. ERROR) GO TO 8 
TEST = 0 .0 

8 CALL ROOTN (TEST ,0.0,2) 

TESTR=P( 1 ) 

T EST I =P ( 2 ) 

DO 18 1=1,2 
S I GN= ( — 1 )**( 1 + 1 ) 

PARTR=EXAM+SIGN*TESTR 

PARTI=SIGN*TESTI 

CALL R00TN(PARTR,PAR7 I ,3) 

IF (TEST .GT. 0.0) GO TO 12 
DO 10 L = 1 ,2 
S I GN= ( - 1 )**(L + 1 ) 

J=2*L-1 
K = 2*L 

AUG ( J ) =P ( 1 ) 

AUG(K)=SIGN*P(2 ) 

10 CONTINUE 
GO TO 20 
12 J = 2 * I - 1 
K = 2*I 

DO 14 L = 1 , 3 

M=2*L-1 

N=2*L 

IP ( ABS ( P ( N ) ) .LE. ERROR) GO TO 16 
14 CONTINUE 
16 AUG ( J ) = P ( M ) 

AUG ( K ) =P ( N ) 

18 CONTINUE 

20 CALL ROOTN (—3 .0,0. 0,2) 

DO 22 1=1,3 
J=2*I-1 
K = 2*I 

ALPHI (J )=AUG( 1 ) +AUG ( 3 ) 

ALPHI (K)=AUG(2)+AUG(4) 



BET I ( J ) =AUG ( 1 ) -AUG ( 3 ) 

BET I ( K ) = AUG ( 2 ) —AUG ( 4 ) 

RG( J ) =P ( 1 ) 

RG( K ) =P ( 2 ) 

22 CONTINUE 
DO 2 A 1=1,3 
J =2* I -1 
K = 2* I 

DEL 1= ( 1-1 ) 

DEL 2= ( 1-2 ) 

DELTA=DEL2/2.0+DELl-6.0*DEL14DEL2/4.0 
DELTB=DEL2/2 .0 

P( J )=-R2/3.0+DELTA*ALPHI ( J ) +DELTB* ( BET I ( J)*RG( JI-BETI (K)*RG(K) ) 
P( K)=DELTA*ALPHI ( K ) +DELTB* ( BET I (K)*RG(J)+BETI(J)*RG(K)) 

24 CONTINUE 
GO TO 38 

26 IF ( ABS ( Z2 ) .LT. 0.00001) GO TO 32 
R 3=Z3/ ( 2 .0*Z2 ) 

R4=Z4/Z2 

TEST=R3**2-R4 

CALL R00TN(TEST,0.0,2) 

DO 28 1=1,2 
J=2*I-1 
K = 2* I 

RG ( J ) =P ( J ) 

R G ( K ) = P ( K ) 

28 CONTINUE 

DO 30 1=1,3 
J=2* I -1 
K = 2* I 

S I GN= ( — 1 ) *# ( 1+1) 

D E L 1 = ( 1-1 ) 

DEL2= ( 1-2 ) 

DELTA=1. O-DEL 1-DEL2/2.0 
P ( J ) =DELTA*(-R3+SIGN*RG ( 1 ) ) 

P(K)=DELTA*SIGN*RG(2.) 

30 CONTINUE 
GO TO 38 

32 IF (Z4.EQ. 0.0) GO TO 36 
DO 34 1=1,3 
J =2* I -1 
K = 2*I 

DEL 1= ( 1-1 ) 

DEL 2= ( 1-2) 

DELT A = 1 .0— DEL l+DELl*DEL2/2.0 
P< J)=-DELTA*(Z4/Z3) 

P ( K ) =0.0 
34 CONTINUE 
GO TO 38 

36 WRITE (6,100) 

WRITE (6,102) 

WRITE (6,104) Z1,Z2,Z3,Z4 

CALL OUT 
CALL EXIT 
38 RETURN 
END 
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i-IBFTC ROOT I LIST, REF, D^CK, DEBUG 

DETERMINATION OF THE ROOTS OF A COMPLEX NUMBER 
SUBROUTINE ROOTN ( A , B , N ) 

COMMON AMO,GAM,THETC,DELB,DELUl,OELU,UC,VC,DELC,BS, DELS, ERROR, 

1 0 , AM , AMU, THETA, A, B,C,D,E,F,G, PRES, DENS, AREA, DYNP, CP, PSI , 

2 XK,YK,UK,VK,XJ,YJ,UJ,VJ,XA,YA,UA,VA,DELA,BETAK,SIGK»SIGSO, 

3 X( 50,50) ,Y( 50,50 ) ,U( 50,50) ,V( 50,50) , BETA! 50) ,PSIR( 50) , 

4 UL,V1,DELL,DELY,RPRES,RDENS,REC0V»ENTP,NDIM,KREF, I REF, 

5 XB( 50) ,YB(50) ,UB( 50) ,VB( 50) ,P( 50) ,RECV(2, 10) ,FACTR, I REV, 

6 REG, R( 25,4 ),S( 25,4 ) ,COWL ( 2 5 ) , BODY ( 25 ) , NR , NS , I COWL, I BODY 
PI=3. 1415927 
EN=N 

EXP= 1 . 0/ EN 

TEST=SORT ( A**2+B**2 ) 

IF (TEST .EO. 0.0) GO TO 16 
AMAG=TEST**EXP 

IF (ABS(A) .GT. 0.0) GO TO 8 
S I GN=B/ ABS ( B ) 

AMAG=ABS(B)**EXP 
AUG= SIGN* PI/2.0 
GO TO 12 

8 IF (ABS(B) .GT. 0.0) GO TO 10 
S I GN = A/ ABS ( A ) 

AMAG=ABS(A)**EXP 
AUG= ( 1 .0— S IGN)*PI/2.0 
GO TO 12 

10 ANGLE = ATAN ( ABS ( B/A ) ) 

SIGNA=A/ABS( A) 

S I GNB = B/ ABS ( B ) 

IF (A .GT. 0.0) GO TO 11 
AUG=P I-S I GNB*ANGLE 
GO TO 12 

11 AUG=( 1.0-SI GNB)*P I +S I GNB* ANGLE 

12 DO 14 1=1, N 
J=2* I — 1 
K = 2* I 

AUG I = ( 1-1 ) 

TAU=EXP*(AUG+2. 0*AUG I *P I ) 

P( J)=AMAG*COS(TAU) 

P(K)=AMAG*SIN(TAU) 

14 CONTINUE 
GO TO 20 
16 DO 18 1 = 1, N 
J=2*I-1 
K=2* I 
P( J)=0.0 
P ( K ) =0 . 0 
18 CONTINUE 
20 RETURN 
END 
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LIST, REF, DECK, DEBUG 


S I BFTC CSAR 
C 

C SUBROUTINE FOR CALCULATING CONIC SHOCK ANGLE 

C 

SUBROUTINE CSA (BETAE ) 

COMMON AMO,GAM, THETC, DEL B,DELU1,DELU,UC, VC, DELC, BS, DELS, ERROR, 

1 Q,AM,AMU,THETA,A,B,C,D, E,F, G,PRES * DENS, AREA ,DYNP,CP, PS I , 

XK, YK,UK , VKtXJt YJ,UJ,VJ,XA, YA, UA, V A , DELA , BET AK , S I GK , S I GSQ , 
X( 50,50) ,Y(50,50) ,U( 50,50) ,V(50,50) ,BETA< 50) ,PSIR(50) , 

U1 , VI t DEL I ,DELY,RPRE$,RDENS,RECOV,ENTP,NDIM,KREF, I REF, 

XB( 50) ,YB<50) ,UB(50) ,VB( 50) t P ( 50 ) , RECV ( 2 , 10 ) , FACTR , I REV , 
REG,R( 25,4) ,S (25,4) ,COWL< 25) , BODY ( 25 ) , NR , NS , I COWL, I BODY 
DIMENSION UU(3) ,VV(3) ,DEL(3) ,ABETA(3) ,TEST(3) ,DBT(3) 

S I GSO=( GAM-1.0) / (GAM+ 1.0) 

00= SORT ( A MO# *2/ { 1 .0+S IGSQ* ( A MO* *2-1 .0 ) ) ) 

DELC=TAN(THETC) 

10 DO 12 1=1,3 
S I GN= 1—2 

ABET A ( I )=BETAE+SIGN*DELB 

Ul=( 1. O-SI GSQ )*QO*COS( ABET A ( I ) ) **2+1. 0/00 
V1 = (00-U1 )*COS(ABETA< I ) )/SIN( ABETAC I ) ) 

DEL 1 = V 1 / { (Jl-QQ ) 

CALL CFF( 1.0, DELC) 

UU( I ) =UA 
VV{ I )=VA 
DEL ( I ) =DEL A 

TEST ( I ) =ABS ( VV ( I )/UU( I ) )-ABS(DELC) 

IF ( ABS ( TEST ( I ) ) .LT. ERROR* ABS ( DELC ) ) GO TO 16 
IF ( ABS ( DELB ) .EQ. 0.0) GO TO 16 
12 CONTINUE 
DO 14 1=1,2 

DBT ( I ) = ( ABETA ( I + 1 J-ABETA ( I ) ) / ( TEST ( 1 + 1 )-TE$T( I ) ) 

14 CONTINUE 

DBDT= ( DBT { 1 ) +D8T ( 2 ) )/2.0 
DELT1=TEST (3) -TEST ( 1 ) 

DELT2=TEST (3)-TEST (2) 

D2BDT2 = ( ( ABETA(3)-ABETA(2) ) -DBDT*DELT2 ) / ( DELT 1 *DELT2 ) 
DELB=-DBDT*TEST ( 2 ) +D2BDT2*TEST ( 1 ) *TEST ( 2 ) 

BETAE=BETAE+ 2.0* DELB 
GO TO 10 

16 0C = SOR T { U1J ( I )**2+VV(I )**2) 

UC=QC*COS ( THETC ) 

VC=OC*SIN(THETC) 

U1=UC 
V1=VC 

DEL 1 =— 1 • 0/ DELC 
BS=ABETA ( I ) 

DEL S = TAN { BS ) 

RETURN 
END 
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SJBFTC CFFR L I ST , REF , DECK , DEBUG 

C 

C FIELD CALCULATIONS IN A CONIC FLOW FIELD 

C 

SUBROUTINE CFF(XP,YP) 

COMMON AMO, G AM, TH ETC » DEL B, DELU 1, OE LU,UC, VC, DELC,BS, DELS, ERROR, 

1 0, AM , AMU, THETA, A, B,C,D,E,F,G, PRES, DENS, AREA, DYNP, CP, PSI, 

? XK,YK,UK,VK,XJ,YJ,UJ»VJ,XA,YA,UA,VA,DELA,BETAK,SIGK,SIGSQ. 

3 X(50,50),Y(50,50) ,U( 50,50) ,V( 50,50), BETA! 50) , PSI R( 50), 

A U1,V1»DEL1»DELY,RPRES»RDENS,REC0V,ENTP,NDIM,KREF, I REF, 

5 XBI50) ,YB(50),UB(50),VB<50) , P ( 50 ) , R EC V ( 2 , 10 ) , F ACTR , I REV , 

6 REG, R < 25,4 ),S( 25,4), COWL (25) , BODY ( 25 ) , NR , NS , ICOWL , I BODY 
DELTA=YP/XP 

T E ST=ABS ( DELTA > - ABS ( 1.0/ DELI ) 

IF (ABS(TEST) .GT. ABS ( DELTA >*ERROR ) GO TO 10 

UA = U1 

VA = V1 

DELA=DEL1 

GO TO 16 

10 SIGNT=TEST/ABS(TEST) 

UP=U1 
VP=V1 
VUP=DEL 1 

DELU=S I GNT*DELU1 

12 CALL RUNGE (UP,VP,VUP,DELU,SIGSO) 

PEST=TEST 

SIGNP=PEST/ABS(PEST) 

TEST=ABS( DELTA) -ABS ( l.O/VUPJ 

IF ( ABS ( TEST ) .LT. ABS ( DELTA ) TERROR ) GO TO 14 
SIGNT=TEST/ABS(TEST) 

IF (SIGNT .EO. SIGNP) GO TO 12 
VUU=FUNV ( UP , VP , VUP , S I GSO ) 

DELtJ = -0 . 5* ( ( 1.0/DELTA )+VUP)/VUU 
IF ( DELU .EO. 0.0) GO TO 14 
GO TO 12 
14 UA=UP 
VA = VP 
DEL A= VUP 
U1=UP 
V1=VP 
DEL 1 =VUP 
16 RETURN 
END 


IBFTC RUNGK LIST, REF, DECK, DEBUG 

RUNGE KUTTA INTEGRATION OF TAYLOR-MACCOLL I EQUATION 

SUBROUTINE RUNGE ( UP , VP , VUP , DU, S I GSQ ) 

REAL K1,K2,K3,K4 
K1=DU*FUNV(UP,VP,VUP,SIGSQ) 

K2=DU*FUNV (UP+DU/ 2.0, VP+DU/ 2.0* VUP+K1*DU/ 8.0, VUP+K 1/8.0, SI GSQ) 
K3=0U*FUNV(UP+DU/2.0,VP+DU/2.0*VUP+Kl*DU/8.0, VUP+K2 / 8 .0, S I GSQ ) 
K4=DU*FUNV (UP+DU,VP+DU*VUP+K3*DU/2.0, VUP+K3 , S I GSQ ) 

UP=UP+DU 

VP=VP+DU*( VUP+(Kl+K2+K3)/6.0) 

VUP=VUP+(Kl+2.0*K2+2.0*K3+K4)/6.0 

RETURN 

END 




$ I BFTC FUN V P LIST t REF T DECK t DEBUG 

C 

C TAYLOR-MACCOLLI EQUATION 

C 

FUNCTION FUNV(UPrVPtVUPfSIGSQ) 
FUN 1 = ( 1.0+VUP**2) 

FUN2 = ( 1 .0-S-IGSQ) s MUP + VP*VUP)**2 
FUN3=< 1.0-SI GSQ*(UP**2+VP**2 ) ) 
FUNV= ( 1 .0/VP)*<FUNl-FUN2/FUN3) 
RETURN 
END 


$ I B FTC CFPRS LISTfREFtDECK t DEBUG 


C 

C 

C 


COMPRESSIBLE FLOW PARAMETER ROUTINE 


SUBROUT 

COMMON 

1 

2 

3 

A 

5 

6 

100 FORMAT 
102 FORMAT 
104 FORMAT 
18HMACH 
106 FORMAT 
108 FORMAT 
110 FORMAT 



Ks i r / / / 

HU, 13X, 1HV, 10X, 


v inu / 

(10X,1P7E14.5//) 

# t i ^ / w / i i w n i nr 


14HVP =1PE12 .5// ) 

V 1 2 FORMAT ( / / 36X ,43HC0MPUTED FLOW FIELD UP TO THE POINT OF EXIT//) 
114 FORMAT ( 1 H 1 ) 

(0= SORT ( CP**2 + VP**2 ) 

AMS0= ( 1. .0-S I GSQ) *0**2 /( 1.0-SI GSO*0**2 ) 

IF ( AMSO .GT. 0.0) GO TO 10 
WRITE (6,102) 

WRITE (6,110) XP,YP,UP,VP 
GO TO 12 
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10 AM= SOR T ( AMS 0 ) 

IF (AM .GT • 1.000) GO TO 16 
WRTTE (6,100) 

WRITE (6,110) XP, YP,UP,VP 
12 CALL OUT 

WRITE (6,114) 

WRITE (6,112) 

WRITE (6,104) 

DO 14 K= 1 , 50 
WRITE (6,106) 

DO 14 1=1,50 

IF ( X ( K , I ) .EQ. 0.0) GO TO 14 
IF (I .GT. K ) GO TO 14 
0=SQRT(U(K,I ) **2 + V ( K , I )**2) 

AM=SORT( (1.0-SIGS0)*Q**2/( 1 . 0-S I GSQ*Q**2 ) ) 

THETA = ATAN(V(K,I ) / U ( K , I ) ) 

WRITE (6,108) X ( K , I ) , Y ( K , I ),0, THETA, U(K, I ) , V ( K , I ),AM 
14 CONTINUE 
CALL EXIT 

16 AMU= ARS I N ( 1 * 0/ AM ) 

THETA= AT AN ( VP/ UP ) 

C = 0/ A M 

A= T AN ( THETA+AMU ) 

B = TAN< THETA-AMU) 

n= -b 

E= -A 

IF (NDIM .GT. 2) GO TO 18 
F = 0 .0 
G = 0 • 0 
GO TO 20 

18 F = 0**2* ( VP/YP )/ ( UP**2-C**2 ) 

G= F 

20 PRES = C**(2.0*GAM/(GAM-1.0) ) 

DENS= C**(2.0/ (GAM-1 .0) ) 

AREA = C**(-2.0/ (GAM-1.0) ) /O 
ENTH= C**2/(GAM-1.0) 

DYNP= DENS*0**2/2.0 

CP= PRES/DYNP 

IF (NDIM .GT. 2) GU TO 22 

PS I =DENS*0/ ( AM*S IN (THETA+AMU) ) 

GO TO 24 

22 PS I =QENS*Q/ 12.0* AM*SIN{ THETA + AMU) ) 

24 RETURN 
END 



o o o 


$ I 8FTC BOUNS LIST » REF, DECK, DEBUG 

CONIC CHARACTERISTIC BOUNDARY ROUTINE 
SUBROUTINE BOUND(KA,IA,KB,IB,XC,YC,SIGR,SIGN) 

COMMON AMO, GAM, THETC, DELB, DELU1,DELU,IJC, VC, DELC,BS, DELS, ERROR, 

1 Q,AM, AMU *THETA,A»B,C,D,E»F,G, PRES, DENS, AREA, DYNP,CP»PSI» 

2 XK,YK,UK,VK,XJ,YJ,UJ,VJ,XA,YA,UA,VA,DELA,BETAK,SIGK,SIGSQ, 

3 X( 50,50) ,Y(50,50) ,U( 50,50) ,V( 50, 50), BETA ( 50) , PSIRI 50) , 

4 U1 , VI ,DEL1 ,DELY,RPRES,RDENS,RECOV,ENTP,NDIM,KREF, I REF, 

5 XB( 50 ) , YB( 50) ,UB( 50) , VB( 50) ,P ( 50 ) , RECV ( 2 , 1 0 ) , FACTR , I REV , 

6 REG,R<25,4) ,S(25,4) ,COWL(25> ,B0DY(25) ,NR , NS , ICOWL , I BODY 
XP=X(KA, I A ) 

YP=SIGN*Y(KA,IA) 

UP=U ( K A , I A ) 

VP=SIGN*V(KA,IA) 

XO=X(KB, IB) 

YQ=SIGN*Y(KB,IB) 

UO=U ( KB , IB ) 

VQ=SIGN*V(KB,IB) 

XR = XC 

YR=S I GN*YC 
SLOPE=SIGN*SIGR 
CALL CFPR(XP,YP,UP,VP) 

AP = A 
DP = D 
F P = F 

CALL CFPR( XQ,YQ»UQ,VQ) 

AREF= . 5* ( A+AP ) 

DREF=0.5*(D+DP) 

FREF=0.5*( F+FP ) 

S I GMA=AREF 
DX=< XQ-XP ) 

DY= ( YO-YP ) 

I F ( ABS ( S I GR ) .GT. ERROR) GO TO 10 
XK = XR 

YK=YP+SIGMA*(XK-XP) 

GO TO 12 

10 XK=(YR-YP+SIGMA*XP—SLOPE*XR)/ (SIGMA- SLOPE) 

YK =YR+ SLOPE* < XK-XR ) 

12 DS=SORT ( DX**2+DY**2 ) 

DV= ( VO-VP ) 

DVDS=DV/DS 

DELSP=SORT ( ( XK-XP)**2+( YK-YP)**2 ) 

VKP=VP+DVDS*DELSP 

UKP=UP+DREF*( VKP-VP )+FREF*( XK-XP ) 

DEL SO=SQRT ( ( XK-XO ) **2+ ( YK-Y0 ) **2 ) 

VKQ= VO— DVDS*OELSQ 

UKO=UQ+DREF*( VKO-VO)+FREF*( XK-XO) 

RATIOP=( DS-DELSP )/0S 
RATI 00= (DS-DELSQ)/DS 
UK=RAT IOP*UKP+RAT I OQ*UKQ 
VK=RATIOP*VKP+RATIQQ*VKQ 
YK = S 1 GN*YK 
VK=S I GN*VK 
RETURN 
END 
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SiBFTC CCSRS L I ST f REF , DECK , DEBUG 

C 

C CONIC CHARACTERISTIC SURFACE ROUTINE 

C 

SUBROUTINE CCSR ( STFR , STF , X I , Y I ,UI t V I , XL »YL t UL f VL) 

COMMON AMO,GAM t THETC,DELB*DELUl,DELU f UC,VC t DELC, BS, DEL 5, ERROR, 

1 Q t AM,AMU f THETA , A ,B,C f D t E f F,G» PRES t D£NS , AREA,DYNP,CP,PSI ♦ 

2 XK t YK,UK,VK ? XJ t YJ f UJ,VJ,XA ,YA f UA,VA t DELA,BETAK f SIGK t $IGSQ, 

3 X(50 f 50) ? Y(50 t 50) tU(50 f 50) f V(50 t 50) t BETA(50) t PSIR(50), 

4 U1 tVI f DELl tDELY t RPRES,ROENStRECOV,ENTP f NDIM,KREF t IREF* 

5 XB<50) t YB<50) ,UB(50) ,VB<50) ,P(50) ,RECV(2 t 10) ,FACTR,IREV, 

6 REG,R(25 t 4) ,S(25,4) f COWL! 25) ,BODY(25) ,NR,NS, ICOWL, I BODY 
CALL CFPR ( X I ,Y I ,UI ,VI ) 

A I = A 
D I =D 
F I = F 

PSII=PSI 

CALL CFPR(XL,YL,UL,VL) 

AL=A 

PS I K = PS I 
A= . S’M A I +AL ) 

10 PSIAV=«5 5 < t ( PSI I + PS IK ) 

IF (NDIM .GT » 2) GO TO 12 
YK=ABS(YI )+(STFR-STF)/PSIAV 
GO TO 14 

12 YK = SORT( YI**2 + < ST FR-ST F ) / PS I AV ) 

14 IF ( Y I .GT. 0.0) GO TO 16 
YK=-YK 

16 XK=XI+( YK-YI )/A 

RATI 0={ YK-YI )/ { Y L — Y I ) 

UK=UI+RATIO*(UL-UI ) 

VK=VI+RATIO*(VL-VI ) 

TEST = PS IK 

CALL CFPRI XK, YK,UK,VK ) 

PS IK=PS I 

IF (ABS(TEST-PSIK) .GT. ERROR ) GO TO 10 

RETURN 

END 


SIBFTC PMERS) L 1ST , REF , DECK , DEBUG 

C 

C PRANDTL MEYER EXPANSION ROUTINE 

C 

SUBROUTINE PMER ( AMP , ANGLE , AMF ) 

COMMON AMO, GAM, THETC, DELB, DELU1, DELU,UC, VC ,DELC,BS, DELS, ERROR, 

1 Q,AM,AMU,THETA,A,B,C,D,E,F,G,PRES,DENS,AREA,DYNP,CP,PSI , 

2 XK,YK,UK,VK,XJ,YJ,UJ,VJ,XA,Ya,UA,VA,DELA,BETAK,SIGK, SIGSQ, 

3 X< 50,50) ,Y( 50,50 ) ,U< 50,50) ,V< 50,50) , BETA <50) ,PSIR<50) , 

4 Ul, VI , DELI, DEL Y , RPRES , ROENS , RECOV,ENTP, NDIM, KREF, I REF, 

5 XB<50) ,YB(50) ,UB(50) ,VB<50) f P ( 50 ) ,RECV ( 2 , 10 ) , FACTR , IREV, 

6 REG,K<25,4) ,S(25,4) , COWL (25) ,BODY(25) ,NR,NS, ICOWL , IBODY 
10 AK= SORT(SIGSO) 

AMFF= SORT ( AMF#*2.0— 1 • 0 ) 

THETAF= { ATAN ( AK*AMFF ) ) / AK-( ATAN { AMFF ) ) 

AMPF= SQRT(AMP**2. 0-1.0) 

THETAP= ( AYAN ( AK*AMPF ) ) / AK-( ATAN < AMPF ) ) 

THETAK= ANGLE + ( THET AP-THETAF ) 

0K= SQRT(AMF**2/(1.0+SIGS0*(AMF**2-1.0) ) ) 

UK=QK* COS(THETAK) 

VK= QK*S I N ( THETAK ) 

RETURN 

END 
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o o o ** oonoo 


LIST, REF, DECK, DEBUG 


$ I BFTC INCS 

INITIAL CHARACTERISTIC ROUTINE 
FOR 2 DIMENSIONS SET DNIM=2 
FOR 3 DIMENSIONS SET DNIM=3 

SUBROUTINE INCR(XI,YI) 

COMMON AMO, GAM, THETC, DEL B, DELU 1 ,DELU,UC, VC, DELC,BS, DELS, ERROR, 

1 0, AM, AMU, THETA, A, B,C,D,E,F,G, PRES, DENS, AREA, DYNP, CP, PSI, 

2 XK,YK,UK,VK,XJ ,YJ,UJ,VJ , X A , Y A , UA , V A , DE L A , B E T AK , S I GK , S I GSO , 

3 X ( 50,50) ,Y< 50,50) ,U( 50,50) , V ( 50, 50 ) , BETA ( 50),PSIR( 50) , 

4 U1,V1,DEL1»DELY,RPRES,RDENS,REC0V,ENTP,NDIM,KREF , I REF, 

5 X B ( 50) ,YB( 50) ,UB( 50) , VB ( 50) ,P( 50) , R EC V ( 2, 10) ,FACTR, IREV, 

6 REG,R(25,4) , S ( 25 ,4 ) ,COWL ( 25 ) ,BODY ( 25 ) , NR , NS , I COWL , I BODY 
T ANT = Y I / X I -DEL Y 
IF (NDIM .GT. 2) GO TO 10 
CALL PSA (AMO, 0.0, THETC ) 

CALL CFPRIXI ,YI ,UK,VK) 

XK=(A#XI-YI ) / ( A — TANT ) 

YK=YI+A*( XK-XI ) 

GO TO 12 

10 DELU=-DELU1 

CALL CFF ( X I , Y I ) 

CALL CFPR( XI ,YI , UA , VA ) 

XK=(A*XI-YI)/(A-TANT) 

YK=YI+A*( XK-XI ) 

CALL CFF ( XK , YK ) 

UK=UA 
VK = VA 

DELU=-DELU1 
12 RETURN 
END 


IBFTC OBSWRS L I S T , R E F , DEC K , DE BUG 
OBLIQUE SHOCK WAVE ROUTINE 
SUBROUTINE OSWR ( AMR, BETAR ) 

COMMON AMO,G AM, THETC, DELB, DELU 1 , DELU, UC , VC, DELC,BS, DELS, ERROR, 

1 0, AM, AMU, THETA, A, B,C,D,E,F,G, PRES, DENS, AREA, DYNP, CP, PSI, 

2 XK,YK,UK,VK,XJ,YJ,UJ,VJ , X A , Y A , UA , VA , DEL A , BET AK , S I GK , S I GSO , 

3 X(50,50) ,Y< 50,50 ) ,U( 50,50) ,V( 50,50 ) ,BETA( 50) , PSI R( 50) , 

4 Ul, VI, DELI, DELY,RPRES,RDENS»RECOV,ENTP, NDIM, KREF,IREF, 

5 XB( 50) , YB( 50) ,UB( 50) ,VB( 50) ,P( 50) ,RECV(2, 10) ,FACTR, IREV, 

6 REG,R ( 2 5,4) , S ( 25 ,4 ) ,COWL ( 25 ) ,BODY( 25) , NR , NS , I COWL , I BODY 
FAMR=( AMR*SIN(BETAR ) )**2 

RPRES=( 2.0*GAM*FAMR-( GAM-1 .0) ) / (GAM+1.0 ) 

RDENS= ( (GAM+1 .0 )*FAMR ) / ( ( GAM-1 . 0 ) *FAMR+2 . 0 ) 

R ECO V = ( R DENS )*#( GAM/ ( GAM— ) .0) )*( 1 .O/RPRES )** ( 1 . 0/ ( GAM-1 .0 ) ) 

IF (RECOV .GT. 0.0) GO TO 10 

ENTP=0.0 

GO TO 5 

10 ENTP = — (GAM-1.0) 4 AL OG (RECOV) 

5 RETURN 
END 
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$ I BFTC PSAS L 1ST, REF, DECK, DEBUG 
PLANE SHOCK ANGLE ROUTINE 
SUBROUTINE PS A ( AMP , THE T P , AL PH A ) 

COMMON AMO, GAM, THETC, DELB, DELU 1 ,DELU,UC, VC, DELC,BS, DELS, ERROR, 

1 0, AM, AMU, THETA, A,B,C, D, E, F , G, PRES , DENS , ARE A, DYNP,CP , PS I , 

2 XK , YK ,UK,VK,XJ,YJ,UJ,VJ,XA,YA,UA,VA,DELA,BETAK,SIGK,SIGSO 

3 XI 50,50) ,Y( 50,50) ,U( 50,50) , V ( 50,50 ) ,BETA( 50) , PS IR< 50) , 

A U1 ,V1 ,DEL1 ,DELY,RPRES,RDENS,RECOV,ENTP,NDIM,KREF, I REF, 

5 XB( 50) ,YB( 50) ,UB(50) ,VB(50) , P ( 50 ) , R EC V ( 2 , 1 0 ) , F AC TR , I RE V * 

6 REG,R(25,4),S(25,4) ,C0WL(25) ,B0DY(25) , NR , NS , I COWL , I BODY 
DIMENSION RG( 10) ,RGI ( 10) 

100 FORMAT ( / / 3 1 X , 62HSUB SON I C CONDITIONS BEHIND THE OBLIQUE SHOCK HAVE 
1 BEEN REACHED// ) 

102 FORMAT ( / / 35X , 53HTHE SHOCK DETATCHMENT CONDITIONS CAN NOT BE SATIS 
1FIED// ) 

104 FORMAT ( //30X,64HTHE SHOCK REFLECTION CONDITIONS AT THE WALL CAN N 
10T BE SATISFIED// ) 

106 FORMAT ( / / 2 5X , 7 1HTHE POINT OF INCIPIENT SHOCK DETATCHMENT HAS BEEN 
1 REACHED, WARNING ONLY//) 

108 FORMAT ( / /45X ,37H I MPROPER INPUT DATA TO SUBROUTINE PSA//) 

110 FORMAT (//46X,3HAMP,10X,5HTHETP,10X,5HALPHA//) 

112 FORMAT ( 40X , 1P3E14*5// ) 

SIGSQ = (GAM-1. 0)/<GAM+1.0) 

AMPSO = AMP*AMP 

0 = SQRT(AMPSQ/(1.0+SIGSQ*(AMPSQ-1.0))) 

IF (ABS(ALPHA) .GT. 0.001) GO TO 8 
UK=0*C0S(THETP) 

VK=0*SIN( THETP ) 

BETAK=ARSIN( 1.0/AMP I+THETP 
S I GK = T AN ( THET P ) 

GO TO 5 

8 CDEFl=-( AMPSQ+2.0)/AMPS0 
C0EF2 = -GAM*SIN< ALPHA ) **2 
C0EF3= ( 2 *0* AMPSQ+ 1.0) / ( AMPSO*AMPSOI 
C OE F4= ( ( GAM+1.0)/2.0)**2+( GAM-1 .0) /AMPSO 
C0EF5=-CGS ( ALPHA ) **2/ ( AMP SO-AMP SO ) 

A A= 1 • 0 

BB=C0EF1+C0EF2 

CC=C0EF3+C0EF4*S IN ( ALPHA )**2 
DD=C0EF5 

10 CALL R00T3I AA,BB,CC,DD) 

DO 12 1=1,3 
J = 2* I 
K = J — 1 

RG( I ) = P ( K ) 

RGI (I )=P( J 
12 CONTINUE 

DO 14 1=1,3 

14 IF ( A8 S ( RG I ( I ) ) . GT * ERROR) GO TO 16 
GO TO 17 

16 WRITE (6,102) 

WRITE (6,104) 

WRITE (6,110) 

WRITE (6,112) AMP, THETP, ALPHA 
CALL OUT 
CALL EXIT 

17 DO 18 1=1,3 

18 IF ( RG ( I ) .LT.O. ) GO TO 20 
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GO TO 21 

20 WRITE (6*108) 

WRITE (6*110) 

WRITE (6*112) AMP*THETP*ALPHA 
CALL OUT 
CALL EXIT 

21 NNR = 0 

22 DO 24 1=1*3 

IF (RG( I ) .GT.O. ) GO TO 24 
NNR=NNR+1 
24 CONTINUE 

IF ( NNR .GT. 0) GO TO 26 
GO TO 28 

26 WRITE (6*108) 

WRITE (6*110) 

WRITE (6,112) AMP*THETP*ALPHA 
CALL OUT 
CALL EXIT 
28 NGR=0 
30 DO 32 1=1*3 

IF( ABS(RG( 1)-RG( I ) >.GT.O. ) GO TO 32 
NGR=NGR+1 
32 CONTINUE 

IF ( NGR • EO • 1 ) GO TO 38 
IF ( NGR • EO. 2 ) GO TO 34 
IF ( NGR • EO *3 ) GO TO 36 
34 WRITE (6*106) 

WRITE (6*110) 

WRITE (6,112) AMP,THETP*ALPHA 
GO TO 38 

36 WRITE (6,108) 

WRITE (6*110) 

WRITE (6,112) AMP*THETP, ALPHA 
CALL OUT 
CALL EXIT 
38 K =0 

HIGH =0 
DO 40 1=1*3 

40 HIGH = AMAX1(HIGH*RG( I ) j 
DO 42 1=1*3 

IF ( RG( I ). EG. HIGH) GO TO 42 
K = K+ 1 

RG(K)=RG( I ) 

42 CONTINUE 
S EC = 0 

DO 44 K = 1 * 2 
44 SEC = AMAX1(SEC*RG(K) ) 

Z = SOR T ( S EC ) 

BET AP = ARS I N ( Z ) 

SINB = SIN(BETAP) 

COSB = COS ( BET AP ) 

U2 = ( (1.0-SIGSQ)*Q*COSB**2)+ 1.0 / 0 
V2 = ( (0 - U2)*C0SB)/$INB 
QK = SORT ( U2**2 + V2**2 ) 

ANGLE=THETP+ALPHA 
UK=QK*COS( ANGLE ) 

VK = QK*SIN( ANGLE ) 

AMK = SORT ( ( ( 1 .0— SIGSQJ^QK# 5 ^ )/ ( l.O-SIGSQ^QK*^ ) ) 
IF (AMK.LE.1.0) GO TO 46 
GO TO 48 
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46 WRITE (6,100) 

WRITE (6,104) 

WRITE (6,110) 

WRITE (6,112) AMP, THETP, ALPHA 
CALL OUl 
CALL EXIT 

48 AMUR = ARS IN ( 1 .0/ AMK ) 
BETAK=BETAP+THETP 
SIGK=TAN( ANGLE) 

5 RETURN 
END 


$ I BETC DRAGS L I ST , REF , DECK , DEBUG 

CALCULATION OF CAPTURE MASS FLOW AND ADDITIVE DRAG IN A CONE FIELD 
SUBROUTINE DR AG ( THET AL , YC ,ND ) 

COMMON AMO, GA M, TH E TC , DEL B, DE LU 1,DELU,UC, VC, DELC,BS, DELS, ERROR, 

1 0, AM, AMU, THETA, A, B,C,D,F,F,G, PRES, DENS, AREA, DYNP, CP, PSI , 

2 XK,YK,UK,VK,XJ,YJ,UJ,VJ,XA,YA,UA,VA,DELA,BETAK,SIGK,SIGSO, 

3 X(50,50) , Y ( 5 0 , 50 ) ,U(50,50) , V ( 50 , 50 ) ,BETA(50) ,PSIR(50) , 

4 U1,V1,DEL1,DELY,RPRES,RDENS,REC0V,ENTP,NDIM,KREF,IREF, 

5 XB ( 50 ) ,YB(50 ) ,IJB( 50) ,VB( 50) ,P< 50) , R EC V ( 2 , 10 ) , F ACTR , I REV , 

6 REG,R (25,4) ,S (25,4) , COWL (25) , BODY ( 25 ) , NR , NS , I COWL , I BODY 
100 FORMAT (//35X,48HINITIAL SHOCK IMPINGES ON INSIDE SURFACE OF COWL/ 

1 / ) 

102 FORMAT ( / / 16X , 1 5HC OWL POSITION = 1 P E 12 . 5 , 2 X , 19HCAPTURE MASS FLOW =1 
1PE12.5,2X,15HADDITIVE DRAG =1PE12.5//) 

104 FORMAT ( / /40X , 39HC0ND I T I ONS ALONG THE CAPTURE STREAMLINE//) 

106 FORMAT ( //8X, 1HX , 1 3X , 1HY , 1 3X , 1HQ, 1 IX , 5HTHET A, 1 1 X, 1 HU, 1 3X , 1HV , 9X, 
18HMACH N0.,8X,4HP/H0,10X,4HP/P0//> 

108 FORMAT (1P9E14.5//) 

110 FORMAT ( 1 H 1 ) 

IF (THETAL .LE. BS) GO TO 10 
WRITE (6,100) 

CALL EXIT 

10 IF (THETAL .LT. BS) GO TO 12 
CAPT=1 .0 
SUM = 0 .0 

WRITE (6,102) BS,CAPT , SUM 
GO TO 42 
12 1=1 

EXP=-(GAM)/ (GAM- 1.0) 

CONV= ( 1 ,0+( GAM-1 .0)/2.0*AMO**2 )**EXP 
CONST = 1 .0/ <GAM*AM0**2 ) 

XND=ND 
T AU = THET AL 
DELT=(BS-TAU)/XND 
XP = YC / T AN ( T AU ) 

YP = YC 
U1=UC 
V1 = VC 

DE L 1 = — 1 .O/DELC 
DELU=DELU1 
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IF (NDIM #GT • 2) GO TO .14 
CALL PSA (AM0,0.0,TF!ETC ) 

UP=UK 
VP = VK 
GO TO 16 

14 CALL C FF ( X P , YP ) 

U P = U A 

V P = VA 

16 CALL CFPR(XP,YP,UP,VP) 

PSI=DENS*Q*SIN(TAU-THETA)/SIN(TAU) 

PRHO=FACTR*PRES 

PRPO = PRHO /CON V 

PS I R ( I ) = PRPO 

WRITE (6,104) 

WRITE (6,106) 

WRITE ( 6 , L08 ) X P , Y P , Q , THET A , UP , VP , AM, PRHO , PR PO 
IF (NDIM .GT. 2) GO TO 18 
STFR=PSI*YP 
GO IU 20 

18 STFR=PSI*YP**2/2.0 
20 SUM = 0 • 0 

DO 36 1=1, ND 
T AU=T AU+DELT 

IF (NDIM . GT . 2) GO TO 22 
UQ=UP 
VO = VP 
GO TO 24 

22 CALL CFF( YC/TANCTAU) ,YC> 

UO=UA 
VQ = VA 

24 CALL CFPR( YC/TAN(TAU) ,YC,UO,VO) 
PSI=DENS*Q*SIN(TAU-THETA)/SIN( TAU) 

IF (NDIM .GT. 2) GO TO 26 
Y Q=ST FR / PS I 
GO TO 28 

26 YO = SORT (2.0*STFR/PSI ) 

28 XO=YO/TAN( TAU) 

CALL CFPR( XO,YQ ,UQ, VO ) 

PRHO=FACTR*PRES 
PR P0 = PRHO/CONV 
PS I R ( 1 + 1 ) =PRPG 

WRITE (6,108) X0,Y0,0,THETA,UQ,V0, AM, PRHO, PRPO 
30 PAVE=(PSIR( I + 1 ) +PS I R ( I ) -2 • 0 ) 

IF (NDIM .GT. 2) GO TO 32 
R EF = YP-Y 0 
GO TO 34 

32 REF=YP**2-YQ**2 
34 PART=CONST*PAVE*REF 
SUM=SUM+PART 

Y P = YO 

36 CONTINUE 

IF (NDIM .GT. 2) GO TO 38 
C APT=YQ 
GO TO 40 
38 CAPT=Y0**2 

40 WRITE (6,102) THETAL ,C APT , SUM 
WRITE (6,110) 

DELS = TAN(THETAL ) 

42 RETURN 
END 


non 


tlBFTC DATUR L I ST , REF , DECK , DEBUG 

STARTING SUBROUTINE FOR OFF DESIGN CALCULATIONS 
SUBROUTINE DATUM ( ST ART , DELO ) 

COMMON AMO, GAM, THETC, DELB, DELU1, DELU,UC» VC, DELC,BS, DELS, ERROR, 

1 0,AM,AMU,THETA,A,B,G..D,E,F,G,PRES,DENS,AREA,DYNP,CP,PSI , 

2 XK,YK,UK,VK,XJ,YJ,UJ,VJ,XA,YA,UA,VA,DELA,BETAK,SIGK,SIGSQ, 

3 X( 50,50) ,Y(50,50) ,U( 50,50) ,V( 50,50) , BETA (50) ,PSIR(50) , 

4 U1,V1,DELI,DELY,RPRES,RDENS,REC0V,ENTP,NDIM,KREF, I REF, 

5 XBI 50) ,YB(50) ,UB(50) ,VB( 50) ,P(50) ,RECV(2, 10) ,FACTR, I REV, 

6 REG, R ( 25,4 ),S< 25,4) ,C0WL(25) ,BODY(25) ,NR,NS, I COWL, I BODY 
100 FORMAT ( //39X,46HINITI AL NET SPACING PARAMETER (START) TO SMALL//) 

K = 1 
1 = 1 

XR=BODY ( 2 ) 

YR=S( 1,1 )+S( 1 ,2>*XR+S ( 1 ,3)#XR**2+S( 1 ,4)*XR**3 
CALL EOUK(K,I,XR,YR,UC,VC> 

SLOPE=A 
UCR=UC 
VCR=VC 
U1=UC 
V1=VC 

DEL1=— 1.0/DELC 
DELU=DELU1 
REF=START/2.0 
K0UT=0 

DO IB K= 1 , 49 

IF (KOUT .EO. 1) GO TO 20 
IF (K .LT. 49) GO TO 10 
WRITE (6,100) 

CALL EXIT 
10 D I ST=START 

IF (K .LE. 2) 01 ST=REE 
ANGLE=ATAN( SLOPE) 

XR=X ( K , I ) +D I ST*COS ( ANGLE ) 

YR=Y(K,I )+DIST*S IN (ANGLE) 

T ANT = YR/ XR 

IF (TANT .EO. DELO) K0UT=1 
IF (TANT .LE. DELO) GO TO 12 
XR=(Y(K,I)-SLOPE*X(K,I ) ) / ( DELO-SLOPE ) 

YR=DELO*XR 
K0UT=1 

12 IF (NDIM .GT. 2) GO TO 14 
UR=UCR 
VR=VCR 
GO TO 16 

14 CALL CFF ( XR, YR ) 

UR = UA 
VR = VA 

16 CALL E0UK(K+1,1,XR,YR,UR,VR) 

SLOPE= A 
18 CONTINUE 
20 RETURN 
END 
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$ I BFTC SPIKEI LIST , REF, DECK , DEBUG 


C 

C 

C 

C 


CHARACTERISTIC ROUTINE FOR ISENTROPIC R AMP , NDIM = 2 
CHARACTERISTIC ROUTINE FOR ISENTROPIC SPIKE* NDIM=3 


SUBROUTINE SP I KE ( K A , I A t M * ST AR T t AMF * S I GN * PR I NT ) 

COMMON AMO,GAM,THETC,DELB,DELUl,DELU,UC,VC,DELC,BS,DELS,ERROR, 


1 

2 

3 

A 

5 

6 

100 FORMAT 
102 FORMAT 
1 0 A FORMAT 


Q , AM , AMU, THET A , A ,B ,C ,0, E , F ,G ,PRES ,DENS , ARE A , DYNP ,CP , PS I , 
XK,YK,UK,VK,XJ, YJ ,UJ , V J , XA , YA,UA, VA, DELA , BETAK , S IGK , S I GSO, 
X( 50,50) ,Y (50,50 ),U( 50,50 ) ,V (50,50) ,BETA( 50) ,PSIR( 50) , 

UL , V 1 , DELI , DELY ,RPRES , RDENS ,RECOV , ENTP ,ND I M, KREF , I REF , 

X6( 50) ,YB(50) ,UB(50) ,VB(50) ,P(50) ,RECV{2,10) ,FACTR, I REV, 
REG, R ( 25, A ),S< 25, A), COWL (25), BODY (25), NR, NS, I COWL, I BODY 
( //A3X,30HC0NDITI0NS ON THE RAMP CONTOUR//) 

( / / A3X ,31HC0NDITI0NS ON THE SPIKE CONTOUR//) 

( //A5X,28HCONDITIONS IN THE FLOW FIELD//) 


106 FORMAT ( / / 8X , 1 HX , 1 3X , 1 HY , 1 3X , 1 HO, 1 1 X , 5HTH E T A , 1 1 X , I HU , 1 3X , 1 HV , 9 X , 
18HMACH NO# ,8X,AHP/H0, 10X ,AHP/P0// ) 

108 FORMAT (1P9E1A.5//) 

110 FORMAT ( 1 H 1 ) 

112 FORMAT (1HJ) 

1 1 A FORMAT ( / / 32 X , 55HC ALCllLAT I ON OF FLOW FIELD PERFORMED IN SUBROUTINE 
1 SPIKE//) 

116 FORMAT (//35X,A6HINITIAL NET SPACING PARAMETER (START) TO SMALL//) 
118 FORMAT (//28X,52HPRESSURE RECOVERY (H/HO) ALONG THE CONTOUR SURFAC 
IE = 1 PE 1 2 • 5// > 


K = K A 
I = I A 


CONVR= ( 1 .0+ (GAM-1 .0) / 2 . 0*AM0**2 ) ** (-GAM/ (GAM- 1.0) ) 
IREV=IREV+1 

R EC V ( 1 , I REV )=RECV ( 1, IREV-1 ) 

RECV(2 , IREV)=RECV<2, I REV-1 ) 

IF (NDIM # GT . 2 ) GO TO 10 
CALL PSA( AM0,0#0,THETC) 

U ( K , I ) =UK 
V ( K , I ) =VK 
UC R = UK 


VCR= VK 
GO TO 12 

10 CALL C FF ( X (K , I ) , Y ( K , I ) ) 

U(K, I ) = U A 
V( K, I )=VA 

12 CALL CFPR( X(K, I ) ,Y(K, I ) ,U(K, I ) ,V(K, I ) ) 
A M P = A M 
ANGLE=THETA 
PS I R ( I )=PSI 
S LOPE = A 
UREF=UC 
VR E F= VC 

DR EF = -1 .O/DELC 
DELU=-DELU1 
K0UT=0 
S TFR= 1 #0 
DO 22 K= 1 , A9 

IF (KOUT .EO# 1 ) GO TO 2A 
IF (K #LT # A9 ) GO TO 1A 
WRITE (6,116) 

CALL CFPR(0.0,0#0,0#0,0.0) 

1A ALPHA=ATAN( SLOPE ) 

XR=X (K,l ) -ST ART ^COS( ALPHA) 
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YR=Y(K, I ) -START*S IN ( ALPHA) 

TANT=YR/XR 

IF (TANT .EO. DELC ) KOUT=l 
IF (TANT .GT. OELC ) GO TO 16 
XR=(Y(K, I )-SLOPE*X(K, I ) ) / ( DE LC-SLOPE ) 

YR=DELC*XR 
KOUT - 1 

16 IF ( NO I M .GT. 2) GO TO 18 
UR=UCR 
VR=VCR 
GO TO 20 
18 U1=UREF 
VI = VREF 
DEL 1 = DR EF 
CALL CFF(XRtYR) 

U R = U A 
VR = VA 
UR EF = U1 
VR E F = V 1 
DREF=DEL 1 

20 CALL E0UK(K+1,I,XR,YR,UR,VR) 

SLOP E = A 

P$IAV=.5*( PSI+PSIR( I ) ) 

EXP = N D I M— 1 

R EF = Y ( K+ 1 , I )**EXP-Y(K, I ) **E X P 
STFR=STFR+PSI AV*REF 
PSIR( I ) =PS I 
22 CONTINUE 
24 KREF=K 
I REF=1 

CALL EOUI (1,X(KREF,1 ) ,Y(KREF,1 ) ,U(KREF,1) ,V(KREF,1 )) 

XM = M 

DEL AM= ( AMP-AMF )/XM 

AMF = AMP 

MR=M+1 

DO 26 1=2, MR 

AMF= AMF- DEL AM 

CALL PMER( AMP, ANGLE, AMF ) 

CALL EOUK ( 1 , I ,X ( 1 , 1 ) , Y ( 1, 1 ) ,UK, VK ) 

CALL CFPR(X(1,I ) , Y ( 1 , I ) , U ( 1 , I ) , V ( 1 , I ) ) 

PS I R ( I ) = PS I 
26 CONTINUE 

DO 38 1 = 1, M 
S TF= 1 .0 

DO 36 K= 1 , K R EF 

CALL CCRC(X(K,I+1 ) , Y { K , I + 1 ) , U ( K , I + 1 ) , V ( K , I + 1 ) , X ( K + l , I ) , Y ( K+ 1 , I ) , 
1U( K+l , I ) , V(K+1 , I ) , S I GN ) 

CALL EOUK (K+l , 1+1 ,XK,YK,UK,VK) 

KRE F = K+ 1 
IREF=I+1 

CALL CFPR(X(K+1,I+1) ,Y(K+1, 1 + 1 ) ,U(K+1,I+1 ) ,V(K+1, I + l) ) 
PSIAV=.5*(PSIR( 1 + 1 )+PSI ) 

IF (NDIM .GT. 2) GO TO 28 
REF=Y(K+ 1,1+1 )-Y(K,I+l) 

GO TO 30 

28 REF=Y (K+l , 1+1 ) **2-Y ( K , I + 1 ) **2 
30 STF=STF+PSIAV*RFF 
PS I R ( 1 + 1 )=PSI 
IF (STF-STFR) 34,32,36 

32 CALL EOU I ( I + l ,X (K+l , 1 + 1 ) ,Y( K + l , 1 + 1 ) ,U(K + 1 , 1 + 1 ) ,V ( K+l , 1 + 1 ) ) 
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GO TO 38 

34 CALL CCSR( STFR ,STF , X ( K + l , 1 + 1 ) , Y ( K+ 1 , I + 1 ) t U ( K + l , I + 1 ) f V ( K + l t 1 + 1 > , 
1X(K,I+1),Y(K,I+1),U(K,I+1),V(K,I+1)) 

CALL EOUI ( I+1,XK.YK,UK,VK) 

GO TO 38 
36 CONTINUE 
38 CONTINUE 

40 CALL EGUK(KREF,IREF,XB( IREF ) ,YB( IREF ), UR ( IREF ) ,VB ( IREF ) ) 

WRITE 16,114) 

IF (NDIM .GT. 2) GO TO 42 
WRITE (6,100) 

WRITE (6,118) RECV ( 1 , I R E V ) 

WRITE (6,106) 

GO TO 44 
42 WRITE (6,102) 

WRITE (6,118) RECV ( 1 , I R EV ) 

WRITE (6,106) 

44 DO 46 1=1 , IREF 

IF ( XB ( I ) .EG* 0.0) GO TO 46 
CALL CFPR ( XB ( I ) , YB ( I ) ,UB( I ) ,VB ( I ) ) 

PRHO=FACTR*PRES 
PRP0=PRH0 /CONVR 

WRITE (6,108) XB( I ) ,YB( I ) ,0, THETA, UB( I ) ,VB( I ) , AM, PRHO , PRPO 
46 CONTINUE 

IF (PRINT .EG. 0.0) GO TO 50 
WRITE (6,110) 

WRITE (6,104) 

WRITE (6,106) 

DO 48 K= 1 , 50 
WRITE (6,112) 

DO 48 1=1, IREF 

IF ( X ( K , I ) .EQ. 0.0) GO TO 48 

CALL CFPR( X(K, I ) ,Y(K, I ) ,U(K, I ) ,V(K,I ) ) 

PRHO=FACTR*PRES 

PRP0=PRH0/C0NVR 

WRITE (6,108) X(K, I ) ,Y(K, I ) ,0,THETA,U(K, I ) , V(K, I) , AM , PRHO , PR PO 
48 CONTINUE 
50 KR = KRE F+ 1 
I R= I R EF — 1 
DO 54 K= 1 , 50 
DO 54 1=1,50 
IF (K .GT. KR) GO TO 52 
IF ( I .GT. IR) GO TO 54 
52 CALL EQUK(K, I ,0.0, 0.0, 0.0, 0.0) 

CALL EOUI ( I, 0.0, 0.0, 0.0, 0.0) 

54 CONTINUE 
KR=KREF+ 1 
DO 56 K =K R , 50 

CALL EQUK (K, IREF, 0.0, 0.0, 0.0, 0.0) 

56 CONTINUE 

DO 58 K= 1 , KREF 
KR=(KREF+1 )-K 

CALL EOUK(K, 1,X(KR, IREF) ,Y(KR, IREF) , U ( KR , I RE F ) , V ( KR , IREF) ) 

CALL EQUK (KR, IREF, 0.0, 0.0, 0.0, 0.0' 

58 CONTINUE 
RETURN 
END 



$ I BFTC CONER LIST , RE F , DECK , DEBUG 

C 

C SUBROUT I Nb FOR CALCULATING THE C-NET BEHIND A SHOCK 

C 

SUBROUTINE CONE ( KA , I A , START , ALPHA , S I GN, PRI NT ) 

COMMON AMO, GAM, THETC, DELB, DELUl,DELU,UC,yC,DELC,BS, DELS , ERROR, 

1 Q, AM, AMU, THETA»A,B,C,D,E,F,G, PRES , DENS ,AR EA ,DYNP ,CP ,PSI , 

2 XK,YK,UK,VK,XJ,YJ,UJ,VJ,XA,YA,UA,VA,DELA,BETAK,SIGK,SIGSO, 

3 X(50,50),Y(50,50),U(50,50),V(50,50),BETA(50),PSIR(50), 

4 U1,V1,DEL1,DELY,RPRES,RDENS,REC0V,ENTP,NDIM,KREF, I REF, 

5 XB ( 50 ) » YB ( 50 ) , UB ( 50 ) ,VB(50) , P ( 50 ) r RECV ( 2 , 10 ) » F ACTR , IREV , 

6 REG,R(25,4),S<25,4) , COWL (25) ,BODY ( 25 ) , NR, NS , ICOWL, I BODY 

100 FORMAT ( //38X,47HC0NDITI0NS IN THE VICINITY OF THE OBLIQUE SHOCK// 
1 ) 

102 FORMAT ( // 1 8X ,8HTAN ( Y/X ) ,6X ,8HMACH NO. , 8X , 5HTHETA ,9X ,4HBETA , 10X , 
13HDEL,10X,5HP2/P1,10X, 5HH2/H1 // ) 

104 FORMAT ( 14X , 1P7E 14. 5// ) 

106 FORMAT ( //43X , 30HC0NDI T I ONS ON THE RAMP CONTOUR//) 

108 FORMAT ( //43X,31HC0NDITI0NS ON THE SPIKE CONTOUR//) 

110 FORMAT ( //45X,2 8HC0NDI T I ONS IN THE FLOW FIELD//! 

112 FORMAT (//8X,1HX,13X,1HY,13X, 1HQ ,11X,5H THETA, 11X, 1HU, 13X,1HV,9X, 
18HMACH NO. ,8X,4HP/H0,10X,4HP/P0// ) 

114 FORMAT ( 1P9E14.5/ / ) 

116 FORMAT ( 1 HI ) 

118 FORMAT ( 1HJ ) 

120 FORMAT ( / /35X, 54HC ALCULAT I ON OF FLOW FIELD PERFORMED IN SUBROUTINE 
1 CONE//) 

122 FORMAT ( //28X , 52HPRESSUKt RECOVERY (H/HO) ALUINU THE CONTOUR SURFAC 
IE =lPE12.5//> 

124 FORMAT (//39X,46H INITIAL NET SPACING PARAMETER (START) TO SMALL//) 
K=KA 
I = I A 

KS I GN=S I GN 

CONVR=( 1.0+ (GAM-1.0 ) /2.0*AM0**2 ) **<— GAM/ (GAM-1. 0) ) 

SUMR=0 .0 
SUMY=0.0 
DELTA=ALPHA 
FACTP=FACTR 

CALL CFPR( 1.0,DELC,UC,VC) 

AMR=AM 

THETAR=THETA 

DELP=TAN(THETAR+ALPHA) 

NREF=NS 

IF (NS .GT. 1) GO TO 10 
X ( K , I )=1.0 
Y ( K , I ) =DELC 

CALL CURVE(0.0,0.0,DELC,X(K,I ) ,0.0 ,0 .0, 1 .0, NS ) 

NS=NS+1 

XP=X ( K , I ) +6.0 

CALL CURVE ( X ( K , I ) , Y ( K , I ) , DEL P, XP, 0. 0, 0. 0 , 1 . 0, NS ) 

NS=NS+1 
10 XS=BODY ( 2 I 

YS=S(2,1 )+S(2,2)*XS+S(2,3) *XS**2+S (2,4)*XS**3 
CALL PSA(AMR,THETC, DELTA) 

CALL EQUK (K,I ,XS,YS ,UK , VK ) 

CALL EOUI ( I ,XS,YS,UK,VK) 

BETA (K )=BETAK 
SLOPE=TAN(BETA(K) ) 

BETAR=BETA(K)-THETAR 
CALL OSWR(AMR,BETAR) 


78 



KREV=(3+KSIGN)/2 

IREV=IREV+1 

RECVtKREV, IREV)=RECV(KREV, I REV-1 )*RFCOV 

R A VE = R EC 0 V 

R EC V 1 = R ECO V 

FACTR=FACTP*RAVE 

WRITE (6,120) 

WRITE (6,122) RECV(KREV, I REV-1 ) 

WRITE (6,100) 

WRITE (6,102) 

WRITE (6,104) DELC , AMR, THE T AR , B E T AR , DELTA, R PRES, R ECOV 
UR E F = UC 
VR EF = VC 

DR EF = - 1 . 0/DELC 
DELU=DELU1 
UC R = UC 
VCR= VC 

KREV=(3-KSIGN)/2 
R E F = S T A RT / 2 • 0 
KOUT = 0 
E XP = ND I M— 1 
DO 32 K= 1 , 49 

IF (KOUT ♦ EO • 1) GO TO 33 
IF (K .IT. 49) GO TO 12 
WRITE (6,124) 

CALL EXIT 
12 DIST=START 

IF (K • L E • 2) D I S T = R E F 
ANGLE=ATAN(SLOPE) 

XR=X ( K , 1 )+DIST*COS(ANGLE ) 

YR=Y(K *1 )+DIST*SIN( ANGLE ) 

TANT=YR/XR 

IF (TANT .EO. DELS) K0UT=1 

IF (TANT .LE. DELS) GO TO 14 

XR= ( Y ( K , 1 )— SLOPE *X(K,1) ) / ( DELS-SLOP E ) 

YR=DELS*XR 

tant=yr/xr 

KOUT = 1 

14 IF ( NO I M .GT. 2) GO TO 16 
UR=UCR 
VR = VC R 
GO TO 18 
16 U 1 = U R E F 
V 1 = VRE F 
DELI =DREF 
DEL U = DEL U 1 
CALL CFF(XR,YR) 

UR = UA 
VR = VA 
UREF=U1 
VR EF = V 1 
DR E F = DEL 1 

18 CALL CFPR ( XR ,YR,UR ,VR ) 

AMR= AM 

THETAR=THETA 

20 IF (K .GT, 1 ) GO TO 22 

CALL CCRA(K,1 ,XR,YR, UR, VR, DELTA, SIGN) 

CALL EQUK ( K+ 1 , 1 ,XK,YK,UK,VK ) 

CALL EOUK (K+1,2,XJ,YJ,UJ,VJ) 

CALL EOUI (2,XJ,YJ,UJ,VJ) 
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BETA (K+l )=BETAK 
SL0PE=TAN(BETA(K+1 ) ) 

BETAR=BETA(K+1 ) -THE TAR 
OELTA = ATAM( V(K+1* 1 )/U(K+l*l > ) -THE TAR 
IF ( ABS ( DELTA ) .LE. ERROR) DELTA=0 • 0 
CALL OS WR ( AMR * 8ETAR ) 

RECV(KREV*IREV)=RECV(KREV*IREV-l)*RECOV 

RBAR=(RECOV+RECVl )/2.0 

YBAR=Y(K+1 * 1 )**EXP-Y (K* 1 )**EXP 

SUMR=SUMR+RBAR*YBAR 

SUMY=SUMY+YBAR 

R AVE = SUMR/ SUMY 

REC Vl-RECOV 

FACTR=FACTP*RAVE 

WRITE (6*104) TANT*AMR*THETAR*BETAR*QELTA*RPRES*RECOV 
CALL C0RE(K+1*2*SIGN) 

GO TO 32 
22 DO 24 J = 2 * 50 

IF ( X ( K * J ) .EQ. 0.0) GO TO 24 

CALL CCRB(K*J*XR*YR*UR*VR* DELTA, SIGN) 

TEST- ( XJ-XK )/ ABS ( XJ-XK) 

IF (TEST .GT. 0.0) GO TO 26 
24 CONTINUE 

26 CALL E0UK(K+1*1,XK*YK*UK*VK) 

CALL E0UK(K+1*J*XJ*YJ*UJ*VJ) 

BETA (K+l ) =BET AK 
SLOPE=TAN( BET A ( K + 1 ) ) 

BETAR = BET A ( K + l ) ~T HET AR 
DELTA=ATAN(V(K+1 *1 )/U(K+l*l) ) -THE TAR 
IF ( ABS ( DELTA ) .LE. ERROR) DEL TA = 0. 0 
CALL OSWR( AMRtBETAR ) 

RECV(KREV t I REV) = RECV ( KREV* I REV-1 ) *RECOV 

RBAR=(REC0V+RECV1 )/2.0 

YBAR=Y(K+1*1 )**EXP-Y(K*1 )**EXP 

SUMR=SUMR+RBAR*YBAR 

SUMY=SUMY+YBAR 

RA VE=SUMR /SUMY 

REC V 1 =RECOV 

FACTR=FACTP*RAVE 

WRITE (6*104) TANT* AMR* THE TAR* BET AR* DELTA *RPRES*RECOV 
DO 28 I = J * 50 

IF ( X ( K * I + 1 ) .EO. 0.0) GO TO 30 

CALL CCRC<X<K*I+1),Y(K*I+1),U(K*I+1)*V(K,I+1)*X(K+1,I)*Y(K+1,I), 
1 U(K+1* I ) * V ( K+l * I ) * -S I GN ) 

CALL EOUK ( K + l * I + 1 * XK * YK *UK * VK ) 

28 CONTINUE 

30 CALL CCRE(X(K+1*I ) , Y ( K+ 1 * I ) *U(K+1*I ) *V(K+1*I ) *SIGN) 

CALL EOUK ( K+l , I+1,XK,YK*UK,VK) 

CALL EOU I ( 1+1 *XK*YK*UK*VK) 

CALL C0RE(K+1* I+l.SIGN) 

32 CONTINUE 

33 K R EF = K 

I REF=KREF 

IF ( NREF .GT. 1 ) GO TO 38 
NS=NREF 

RATI0=1.0/Y (KREF* 1 ) 

DO 34 1=1*50 

CALL EOUK I * RAT I G*XB ( I ) * RATIO* YB( I ) *UB( I ) * VB { I ) ) 

DO 34 K = 1 * 50 

IF ( X ( K * I ) .EO. 0.0) GO TO 34 
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CALL EQUK(K,I*RATIO*X(K»I ) , RAT IO*Y ( K , I > ,U(K, I) ,V(K, I ) > 

34 CONTINUE 

CALL CURVE(0.0,0.0,DELC,X< 1 ,1), 0.0,0. 0,1.0, NS) 

NS=NS+1 

XP=X(KREF,IREF)+3.0 

36 CALL CURVE(X(l,l),Y(l,l),V(l,l)/U(l,l),XP,0.0,0.0,1.0,NS) 
NS=NS+1 

38 WRITE (6,116) 

40 IF (NDIM .GT. 2) GO TO 42 
WRITE (6,106) 

WRITE (6,112) 

GO TO 44 
42 WRITE (6,108) 

KREV= ( 3+KS IGN ) t'c 

WRITE (6,122) RECV ( KREV , I REV ) 

WRITE (6,112) 

44 DO 48 1=1,50 

IF ( XB ( I ) .EO. 0.0) GU IU 48 
46 CALL CFPR(XB( I ) ,YB( I ) ,OB( I ) ,VB( I ) ) 

PRHO=FACTR*PRES 

PRPO=PRHO/CONVR 

WRITE (6,114) XB( I ) ,YB( I ) ,0, THETA, UB( I ) ,VB( I ) ,AM, PRHO, PRPO 
48 CONTINUE 

IF (PRINT .EO. 0.0) GO TO 54 
WRITE (6,116) 

WRITE (6,110) 

WRITE (6,112) 

50 DO 52 K=1 , 50 
WRITE (6,118) 

DO 52 1=1,50 

IF ( X ( K , I ) .EO. 0.0) GO TO 52 
IF ( I .GT. K) GO TO 52 

CALL CFPR(X(K,I ) , Y ( K » I ) *U( K , I ) , V ( K , I ) ) 

PRHO=FACTR*PRES 

PRPO=PRHO/CONVR 

WRITE (6,114) X(K,I),Y(K,I ) , 0, THET A , U( K , I ) , V ( K , I ) , AM, PRHO, PRPO 
52 CONTINUE 
54 RETURN 
END 


non 


$ I BFTC PUTS LIST, REF, DECK, DEBUG 

SUBROUTINE FOR LOCATING THE COWL CONDITIONS 
SUBROUTINE PUT ( XCOWL , YCOWL ) 

COMMON AMO, GAM ,THETC, DELB, DELU1 ,DELU,UC, VC, DELC,BS, DELS, ERROR , 

1 0,AM,AMU,THETA,A,B,C,0,E,F,G,PRES,DENS,AREA,DYNP,CP,PST t 

2 XK,YK,UK,VK,XJ,YJ,UJ',VJ,XA,YA,UA,VA,DELA,BETAK,SIGK,SIGSO, 

3 X( 50,50) ,Y( 50,50) ,U( 50,50) ,V( 50,50) ,BETA( 50) ,PSIR( 50) , 

A U 1 , V 1 , DELI, DELY,RPRES,R DENS, RECOV,ENTP,NDIM,KREF,IREF, 

5 XBI50) ,YB( 50) ,UB( 50) ,VB( 50) ,P( 50) ,RECV(2, 10) ,FACTR, I REV, 

6 REG,R ( 25,4 ) ,S( 25,4) , COWL (25) , BODY (25) , NR , NS , I COWL , I BODY 
100 FORMAT ( //25X,57HC0WL LIP POINT LIES OUTSIDE REGION OF FLOW FIELD, 

1 XCOWL =1PE12.5,2X,7HYC0WL =1PE12.5//) 

DO 10 K=1 , 50 

IF ( X ( K , 1 > .EO. 0.0) GO TO 12 
10 CONTINUE 
12 KP=K— 1 

I) I ST = Y ( KP, 1 ) — YCOWL 
IF (DIST .GE. 0.0) GO TO 14 
WRITE (6,100) XCOWL, YCOWL 
CALL CFPR(0. 0,0. 0,0. 0,0.0) 

14 DO 16 1=1, KP 

IF ( X ( KP , I ) .EO. 0.0) GO TO 16 
DIST=X(KP,I )- XCOWL 
IF (DIST .GE. 0.0) GO TO 18 
16 CONTINUE 

WRITE (6,100) XCOWL, YCOWL 
CALL CFPR(0. 0,0. 0,0. 0,0.0) 

18 RM I N= 10 .0 
DO 20 K= 1 , 50 
DO 20 1=1,50 

IF ( X ( K , I ) .EO. 0.0) GO TO 20 
DELX=X ( K , I ) -XCOWL 
DEL Y=Y ( K , I )— YCOWL 
RADIUS=S0RT(DELX**2+DELY**2) 

RMIN=AMIN1 (RADIUS,RMIN) 

20 CONTINUE 

DO 22 K=1 , 50 
DO 22 1=1,50 

IF ( X ( K , I ) .EO. 0.0) GO TO 22 
DELX=X ( K , I )— XCOWL 
DEL Y=Y ( K , I ) -YCOWL 
RADIUS=S0RT(DELX**2+0ELY**2) 

IF (RADIUS .EO. RMIN) GO TO 24 
22 CONTINUE 

24 IF ( ABS ( RADIUS ) .GT. ERROR) GO TO 26 
KREF=KP 
I REF= I P 
GO TO 38 
26 KP=K 
I P= I 

KREF=K+1 

CALL CFPR(X(KP,IP),Y(KP,IP),U(KP,IP),V(KP,IP)) 

TANA=A 
T ANB = B 
DO 28 J=1 ,2 
KSIGN=<-1)**< J+l ) 

SIGN=(-1 )**( J+l ) 

KREP=KP+KS I GN 
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IF (X(KREPflP) * EQ* 0.0) GO TO 28 

CALL BOUND(KP, IP,KREP, IP,XCUWL, YCOWL, TANB, 1.0) 

TEST=(XK-X(KP, IP) )/ABS(XK-X(KP, IP) ) 

IF (TEST .EQ. SIGN) GO TO 30 
28 CONTINUE 

30 X I =XK 

Y I = YK 
U I =UK 

V I = VK 

DO 34 J = 1 t 2 
I R E P= I P 

ISIGN*(-1 )**( J+l ) 

IF (IP .EQ. 1) I S I GN= 1 
S I GN= I S IGN 

31 IREP=IREP+ISIGN 

IF (X(KPflREP) .GT. 0.0) GO TO 32 
GO TO 31 

32 CALL B01JND(KP, IP, KP t IREP,XCOWL, YCOWL, TANA, -1.0) 
1EST=( XK-X(KP, IP) )/ABS(XK-X(KP,IP) ) 

IF (TEST .EQ. SIGN) GO TO 36 
34 CONTINUE 

36 CALL CCRC(XK,YK,IJK,VK,XI,YI,UI,VI,1.0) 
kref=kp 

I REF= IP 

CALL EQUK(KREF f I REF, XCOWL , YCOWL , UK , VK ) 

38 RETURN 
END 


SIBFTC SHAPS L I ST , REF , DECK t DEBUG 

C 

C CHARACTERISTIC FIELD FOR THE BOUNDARY Y=A0+A1* Z + A2* Z**2 + A'3*Z**3 

C 

SUBROUTINE SHAPE ( K A , I A , S I GN , S ET , PR I NT ) 

COMMON A MO, GAM, THETC , DEL B , DE LU 1 , DELU , UC , VC , DELC , BS , DEL St ERROR, 

1 Q,AM,AMU,THETA,A,B,C,D,E,F,G,PRES,DENS,AREA,DYNP,CP,P$I f 

2 XK, YK,UK, VK,XJ, YJ,UJ, VJ,XA, YA,UA, VA,DELA,8ETAK, SIGK,SIGSQ, 

3 X( 50,50) ,Y(50, 50) ,U( 50,50) ,V(50,50),BETA( 50) ,PSIR (50), 

4 U1 ,V1,DEL1,DELY,RPRES,RDENS,REC0V,ENTP,NDIM,KREF, I REF, 

5 XB(50) ,YB(50),UB(50) ,VB(50) , P ( 50 ) , R ECV ( 2 , 10 ) , F ACTR , I REV , 

6 REG, R( 2 5,4 ) ,S( 2 5,4) , COWL ( 25) , BODY (25) ,NR,NS, I COWL, I BODY 

100 FORMAT ( //30X.30HC0ALESCENCE HAS OCCURED AT X =1 PEI 2 . 5 , 2 X , 3HY = 1PE 
112.5// ) 

102 FORMAT ( / /43X , 33HC0NDI T IONS ON THE CONTOUR SURFACE//) 

104 FORMAT ( / /45X , 28HCQNDI T IONS IN THE FLOW FIELD//) 

106 FORMAT ( / / 8X , 1 HX , 1 3X , 1HY , 1 3 X , 1HQ , 1 1 X , 5HTHETA , 1 1 X t 1HU , 1 3X , 1HV , 9X , 
18HMACH NO.,8X,4HP/HO,10X,4HP/PO// ) 

108 FORMAT (1P9E14.5//) 

110 FORMAT ( 1 HI ) 

112 FORMAT (1HJ) 

114 FORMAT ( // 33X , 55HC ALCUL AT I ON OF FLOW FIELD PERFORMED IN SUBROUTINE 
1 SHAPE//) 
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116 FORMAT (//28X,52HPRESSURE RECOVERY (H/HO) ALONG THE CONTOUR SURFAC 
IE =1PE12.5//) 

WRITE (6,110) 

WRITE (6,114) 

CONVR=( 1.0+ (GAM-1 . 0 ) /2 . 0*AM0**2 ) **( -GAM/ (GAM-1.0) ) 

KS I GN=S I GN 
I R E V= I R E V+ 1 

REC V ( 1 , I REV ) = R£CV ( 1 , I REV-1 ) 

RECV(2, IREV)~RECV(2, IREV-1 ) 

12 K = K A 
I = I A 

CALL EOUI'( I ,X(K, I ) ,Y(K, I ) ,U(K, I ) ,V(K, I ) ) 

KREF=KA 
I R= I A+ 1 

00 30 I = I R , 40 
KREF=KREF+1 

K = KREF 
LP=I-1 

DO 14 L = 1 , L P 
J = I-L 

IF ( X ( K , J ) .GT. 0,0) GO TO 16 
14 CONTINUE 
GO TO 32 

16 CALL CCRE(X(K,J) ,Y(K,J) ,U(K,J) ,V(K,J),SIGN) 

CALL EOUK (K,K,XK,YK,UK,VK) 

CALL EOUI (K,XK,YK,UK,VK) 

1 REF= I 

DO 28 K=K R EF , 50 
J R= I- I A 
DO 18 J = 1 , JR 
I PT = I— J 

IF ( X ( K+ 1 , I PT ) .GT. 0.0) GO TO 20 
18 CONTINUE 
GO TO 30 

20 TEST=X(K+1, I PT ) — X ( K , IPT) 

IF (ABS(TEST) .GT. 0.0) GO TO 22 

CALL EOUK ( K+ 1 , I,X(K,I),Y(K,I),U(K,1KV(K,I)) 

GU I U 24 

22 CALL CCRC(X(K+1,IPT),YC<+1,IPT),U(K+1,IPT),V(K+1,IPT), 

1 X(K,I),Y(K,I),U(K,I),V(K,I),SIGN) 

CALL EOUK ( K+l , I ,XK,YK,UK,VK) 

RADII=X(K+1, I ) -X ( K , I ) 

IF (RADI I .GT. 0.0) GO TO 24 
WRITE (6,100) X(K+1 , I) ,Y(K+1, I ) 

IF (K .EO. I ) GO TO 30 

CALL EOUK (K, I ,X (K+l, I ) , Y ( K+l , I ) , U ( K+ 1 , I) , V ( K+ 1 , I ) ) 

24 RADIU$=X(K+1, I )-X(K+l,IPT) 

IF (RADIUS .GT. 0.0) GO TO 28 
WRITE (6,100) X ( K+l , I ) , Y ( K+l , I ) 

IF (IPT .EG. IA) GO TO 30 
KF=K+ 1 

DO 26 J PT = K F , 50 

CALL EOUK ( J PT , IPT, 0.0, 0.0, 0.0, 0.0) 

26 CONTINUE 
28 CONTINUE 
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30 CONTINUE 
32 WRITE (6,102) 

KR E V = ( 3+KS IGN) /2 

WRITE (6,116) RECV ( KRE V , I RE V ) 

WRITE (6,106) 

DO 34 1=1,50 

IE (XB( I) • E.Q • 0.0) GO TO 34 
CALL CFPR(XB( I ) ,Y6( I ) ,UB( I ) ,VB< I ) ) 

PRHO=FACTR*PRES 

PRPO=PRHO/CONVR 

WRITE (6,108) XBU ) 9 YBU > ,0,THETA,UB(I ) ,VB(I ) , AM , PRHO , PR PO 
34 CONTINUE 

IF (PRINT .EO. 0.0) GO TO 37 
WRITE (6,110) 

WRITE (6,104) 

WRITE (6,106) 

DO 36 K=1 ,KREF 
WRITE (6,112) 

DO 36 1=1,50 

IF { X ( K , I ) .EO. 0.0) GO TO 36 
IF ( I .GT. K ) GO TO 36 

CALL CFPR(X(K,I ),Y(K,I ),U(K, I ),V(K,I ) ) 

PRHO=FACTR*PRES 

PRPO=PRHO/CONVR 

WRITE (6,108) X(K,I ),Y(K,I) , 0 , THETA , U ( K , I ) , V ( K , I ) , AM , PRHO , PR PO 

36 CONTINUE 

37 DO 38 K =K A , 50 

IF ( X ( K , I A ) .EQ. 0.0) GO TO 40 

38 CONTINUE 
40 KREF=K-1 

IF (SET .EO. 0.0) GO TO 56 
K P = KREF 
IP = 1 

CALL MATRX 

CALL ENDS(KP,1 SIGN) 

KRE F = K P 
I R E F= I P 

CALL E01JK(KREF,KREF,XB (KREF ) , YB ( KR E F ) , UB ( KR E F ) , VB ( K R EF ) ) 

K = 0 

DO 46 J = 1 , 50 

IF ( X ( J ,KR E F ) .EO. 0.0) GO TO 46 
K = K+ 1 

CALL E0UK(K,1 ,X( J ,KREh) , Y ( J , KR EF ) , U ( J , KR EF ) , V ( J , KRE F ) ) 

46 CONTINUE 
48 KR EF = K 

DO 50 K= 1 , 50 
DO 50 1=2,50 

CALL EOUK(K, I ,0.0,0.0,0.u,0.0) 

50 CONTINUE 

CALL E0UI(1,X(1,1),Y(1,1),U(1,1),V(1,1) ) 

DO 52 1=2,50 

CALL EOUI ( I ,0.0,0. 0,0. 0,0.0) 

52 CONTINUE 
56 RETURN 
END 
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$ I BFTC SKOCKS L 1ST, REF, DECK, DEBUG 

SU8ROUT t NE FOR CALCULATING THE C-NET BEHIND A SHOCK 

FOR INTERNAL ISEN. COMP. ON DESIGN SETP=0.0, SET0=0.0 
FOR INTERNAL ISEN. COMP. OFF DESIGN SETP=1 ,0,SET0=1 .0 
FOR INTERNAL REFLCT . SHOCK ON DESIGN SETP=1.0,SETQ=0.0 
FOR INTERNAL REFLCT. SHOCK OFF DESIGN SETP=1 .0,SET0=0.0 

SUBROUTINE SHOCK (KA, l A, DELTA, SIGN, SETP,SETO) 

COMMON AM0,GAM,THETC,DELB,DELU1,DELU,UC,VC,DELC,BS, DELS, ERROR, 

1 0, AM , AMU, THETA, A, B,C,D,E,F,G, PRES, DENS, AREA, DYNP, CP, PSI, 

2 XK,VK,UK,VK,XJ,YJ,UJ,VJ,XA,YA,UA,VA,DELA,BETAK,SIGK,SIGSO, 

3 X( 50,50) ,Y(50,50) ,U< 50,50) ,V( 50,50) ,BETA( 50) ,PSIR( 50) , 

4 U1,VL,DEL1,DELY,RPRES,RDENS,REC0V,ENTP,NDIM,KREF, I REF, 

5 XB( 50) »YB(50) »UB(50) ,VB( 50) , P ( 50 ) , R EC V ( 2 , 10 ) , FACTR , I REV , 

6 REG,R(25»4) ,S(25,4) ,C0WL(25) , BODY (25), NR, NS, ICOWL, I BODY 
100 FORMAT ( //38X,47HC0NDITI0NS IN THE VICINITY OF THE OBLIQUE SHOCK// 

1 ) 

102 FORMAT ( / / 1 5X , 1HX , 1 2X , 1HY , 10X , 8HMACH NO . , 8X , 5HTHETA ,9X , 4HBETA , 1 OX , 
13HDEL, 10X,5HP2/P1, 10X,5HH2/H1//) 

104 FORMAT ( 7X , 1 P8E 14.5/ / ) 

106 FORMAT ( / / 8X , 1HX , 1 3X , 1HY , 1 3X , 1H0, 1 1 X , 5HTHETA, 1 1 X , 1 HU, 1 3X , 1 HV , 9X , 
18HMACH NO. ,8X,4HP/H0» 10X ,4HP / PO/ / ) 

108 FORMAT (1P9E14.5//) 

110 FORMAT (//43X,33HC0NDITI ONS ON THE CONTOUR SURFACE//) 

112 FORMAT (//45X,28HC0NDITI ONS IN THE FLOW FIELD//) 

114 FORMAT ( 1 HI ) 

116 FORMAT (1HJ) 

118 FORMAT ( //34X,55HCALCULATI0N OF FLOW FIELD PERFORMED IN SUBROUTINE 
1 SHOCK//) 

120 FORMAT ( //28X,52HPRESSURE RECOVERY (H/HO) ALONG THE CONTOUR SURFAC 
IE =1PE12.5//) 

WRITE (6,114) 

WRITE (6,118) 

WRITE (6,100) 

WRITE (6,102) 

CONVR= ( 1 .0+ ( GAM-1 .0 ) /2 ,0*AM0**2 ) *# ( -GAM/ ( GAM-1.0) ) 

KS I GN=S I GN 
SUMR=0.0 
SUMY=0 . 0 
DEL=DELTA 
FACTP=FACTR 

CALL CFPR(X( IA,KA) ,Y( IA,KA) ,U< I A ,KA ) , V ( I A , KA ) ) 

CALL PSA(AM, THETA, DEL) 

CALL EQUK( 1 ,1 ,X( I A ,K A ) , Y ( I A ,KA ) ,UK , VK ) 

CALL E QU I ( 1 » X ( I A , K A ) , Y ( I A , K A ) , UK , V K ) 

CALL CFPR(X(IA,KA),Y(IA,KA),U(IA,KA),V(IA,KA)) 

AMR=AM 

BETA ( 1 )=-SIGN*BETAK+< 1 ,0+S I GN ) *THETA 
SLOPE = TAN(BETA( 1 ) ) 

BETAR=BETA(1)-THETA 
KREF=KA 
I R EF= I A 

CALL OSWR(AMR,BETAR) 

KREV=(3+KSIGN)/2 
I RE V= I R E V+ 1 

RECV(KREV, IREV)=RECV(KREV, I REV-1 )*RECOV 
RAVE=RECOV 
RECV1=REC0V 
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FACTR=FACTP*RAVE 

WRITE (6,104) X( IREF, KREF), Y( IREF»KREF) , AMR , THETA , BET AR ,DEL , R PRES, 
1REC0V 

KREV=(3-KSIGN>/2 
EXP=ND I M— 1 
DO 30 K=2 , 50 

IF (KREF ,E0. 50) GO TO 32 

IF ( ABS ( DELTA ) .GT. 0.001) GO TO 10 

1 REF= I R EF+1 

IF ( X ( IREF ,KREF ) .EO. 0.0) GO TO 32 
GO TO 14 
10 DO 12 L= 1 , 2 
REG=— REG 

CALL SEEK (KREF, I REF, X (I REF, KREF ) ,Y( IREF, KREF) , SLOPE, SIGN) 

CALL EOUK ( IREF, KREF, XK,YK, UK, VK) 

REG=-REG 

IF (K .EO. 2) GO TO 12 
IF (KREF .EO. 50) GO TO 14 
DEL X=X ( K-l , 1 ) -X ( I REF , KREF ) 

DELY=Y(K-1,1 ) — Y ( IREF, KREF) 

DIST=SORT(DELX**2+DELY**2) 

IF (ABS(DIST) .GT. 0.000) GO TO 14 
12 CONTINUE 

14 IF (K .GT. 2) GO TO 20 

16 CALL CCRA < K— 1 , 1 , X ( I REF ,KREF ) , Y ( I REF , KREF ) , U ( I R EF , KREF ) , 

1 V( IREF, KREF) , DEL, SIGN) 

CALL EOUK (K,1,XK»YK,UK,VK) 

CALL EOUK (K»2,XJ,YJ»UJ»VJ) 

CALL EOUI (2»XJ,YJ,UJ,VJ) 

CALL C FPR ( X ( IREF, KREF) ,Y ( IREF, KREF) ,U ( I R EF ,KRE F ) , V ( IREF, KREF) ) 
AMR=AM 

BETA (K)=-SIGN#BETAK+( 1.0+ SIGN )*THETA 
SLOPE = TAN(BETA(K ) ) 

BETAR=BETA(K)— THETA 
DEL=ATAN(V(K»1)/U(K,1) (-THETA 
IF ( ABS ( DEL ) .LE. ERROR) DEL=0.0 
CALL OSWR(AMR,BETAR) 

RECV(KREV, IREV)=RECV(KREV, I REV-1 )*RECOV 
RBAR= ( REC0V+RECV1 ) /2 .0 
YBAR=ABS(Y(K,1 ) **EXP-Y (K-l , 1 ) **fXP ) 

SUMR=SUMR+RB AR*YB AR 

SUMY=SUMY+YBAR 

RAVE=SUMR/SUMY 

RECV1=REC0V 

FACTR=FACTP*RAVE 

WRITE (6,104) X( IREF, KREF),Y( IREF, KREF) , AMR , THETA , BET AR , DEL , R PR ES , 
1REC0V 

CALL CORE ( K ,2 , S I GN ) 

GO TO 30 
20 DO 22 J =2 , 50 

IF ( X ( K-l , J ) .EO. 0.0) GO TO 22 

CALL CCRB ( K — 1 , J , X ( I REF , KREF ) , Y ( I REF ,KREF ) ,U t I REF ,KREF ) , 

1 V( IREF, KREF) , DEL, SIGN) 

TEST=(XJ-XK )/ABS(XJ-XK) 

IF (TEST .GT. 0.0) GO TO 24 
22 CONTINUE 

24 CALL EOUK ( K , 1 , XK , YK , UK , VK ) 

CALL EOUK(K,J,XJ,YJ,UJ,VJ) 

K P=K 

CALL C FPR ( X ( IREF, KREF ) » Y ( I REF, KREF) ,U{ IREF, KREF), V( IREF, KREF)) 
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AMR= AM 

BETA (K ) =-S IGN*BETAK+( 1 . 0+ S I GN ) *THE T A 
SLOPE=TAN(BETA(K) ) 

BETAR=BETA (K )-THETA 
DEL=ATAN(V(K»1)/U(K,1) )-THETA 
IF (ABS(OEL) .LE. ERROR) DEL=0.0 
CALL OSWR ( AMR,BETAR ) 

RECV(KREV,IREV)=RECV(KREV,IREV-1 )*RECOV 

RBAR=(REC0V+RECV1 )/2.0 

YBAR=ABS( Y (K ,1 ) **EXP-Y( K-l , 1 ) #*EXP > 

SUMR=SUMR+RBAR*YBAR 
SUMY=SUMY+YBAR 
RAVE=SU MR/SUMY 
RECV1=REC0V 
FACTR=FACTP*RAVE 

WRITE (6,104) X ( IREF, KREF ) , Y ( I REF ,KREF ) , AMR , THETA, BET AR,DEL,R PRES, 
1REC0V 

00 26 I = J , 50 

IF ( (K-l ) .EO. I ) GO TO 28 

CALL CCRC(X(K,I ) ,Y(K, I ) ,U(K,I ) ,V(K,I ) ,X(K — 1,1 + 1) ,Y(K-1, 1+1 ) , 

1 U(K-1 , 1+1 ) , V (K-l, I + 1 ) , S I GN ) 

CALL E0UK(K,I+1,XK,YK,UK,VK) 

26 CONTINUE 

28 CALL CCRE ( X (K , I ) ,Y ( K, I ) ,U( K, I ) ,V (K , I ) , S I GN ) 

CALL EOUK(K, I+1,XK,YK,UK,VK) 

CALL EOUI ( 1 + 1 ,XK,YK,UK,VK) 

CALL CORE(K, I+1,SIGN) 

IP=I+1 
30 CONTINUE 
32 KREF=KP 
I REF= I P 

IF ( SETP .GT. 0.0) GO TO 34 
CALL LOCUS ( 1 , S I GN ) 

STFR=PSI 

34 WRITE (6,114) 

WRITE (6,110) 

KREV=(3+KSIGN)/2 

WRITE (6,120) RECV ( KREV , I REV ) 

WRITE (6,106) 

00 36 1=1, IREF 

IF ( XB ( I ) .EO. 0.0) GO TO 36 
CALL CFPR(XB( I ) ,YB( I ) ,UB( I ) ,VB( I ) ) 

PRHO=FACTR*PRES 

PRPO=PRHO/CONVR 

WRITE (6,108) X6( I ) ,YB ( I ) ,0, THETA, UB( I ) ,VB( I ) ,AM,PRHO,PR 
36 CONTINUE 

WRITE (6,114) 

WRITE (6,112) 

WRITE (6,106) 

DO 38 K=1 ,KREF 
WRITE (6,116) 

DO 38 1=1,50 

IF ( X ( K , I ) .EQ. 0.0) GO TO 38 
IF (I .GT. K) GO TO 38 

CALL CFPR(X(K, I ),Y(K,I ),U(K, I ),V(K,I ) ) 

PRHO=FACTR*PRES 

PRPO=PRHO/CONVR 

WRITE (6,108) X(K, I ) ,Y(K, I ) ,0 , THETA ,U ( K , I ) ,V(K,I ) ,AM,PRHO,PRPO 
38 CONTINUE 
KF=KREF+1 
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00 46 K=KF , 50 
00 46 1=1,50 

CALL EOUMK, I ,0.0,0.0,0.0,0.0) 

46 CONTINUE 

DO 50 1=1,50 

IF (I .GT. KREF) GO TO 48 

CALL EOUI ( I ,X( I , I ),Y< I , I ) ,U< 1,1) , VII , I ) ) 

GO TO 50 

48 CALL EOUI ( I ,0.0,0. 0,0. 0,0.0) 

50 CONTINUE 

IF ( SETP .EO. 0,0) GO TO 62 
KP=KREF 
I P= I R E F 
CALL MATRX 

CALL ENDS ( KP , 1 *-S I GN ) 

KREF=KP 
I R EF= I P 

CALL EOUK ( KREF, KREF, XB( KREF) ,YB< KREF) ,UB< KREF) ,VB( KREF) ) 

IF (SETO . EO . 0.0) GO TO 62 

K=0 

DO 54 J = 1 , 50 

IF (X(J,KREF) ,E0. 0.0) GO TO 54 
IF (J .EO. 1) GO TO 52 
DELX=Xt J,KREF)-X(K,1) 

DELY=Y(J,KREF)-Y(K,1) 

DIST=SQRT(DELX**2+DELY**2 ) 

IF (ABS(DIST) .LT. 0.000) GO TO 54 
52 K =K+ 1 

CALL EOUK ( K, 1, X ( J, KREF ),Y(J, KREF) , U ( J , KREF ), V ( J , KREF ) ) 

54 CONTINUE 
KF = K 

DO 58 K=1 , 50 
DO 58 1=2,50 

56 CALL EOUK (K,I ,0.0,0. 0,0. 0,0.0) 

58 CONTINUE 

CALL EOUI (1,X(1,1),Y(1,1),U(1,1),V(1,1)) 

DO 60 1=2,50 

CALL EOUI (I ,0.0, 0.0, 0.0, 0.0) 

60 CONTINUE 
62 PSI=STFR 
DEL A=DEL 
RETURN 
END 
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$1BFTC SEEKS L I ST , RE F T DECK , DEBUG 

SUBROUTINE FOR LOCATING SPECIFIC POINTS IN A C-NET 
SUBROUTINE SEEK ( K A , I A , XR , YR , SLOPE , S I GN ) 

COMMON AMO, GAM, THETC, DELB, DELU1,DELU,UC, VC, DELC, 8 S, DELS, ERROR, 

1 0, AM, AMU, THETA, A, B,C,D,E,F,G, PRES, DENS, AREA, DYNP, CP, PSI , 

2 XK, YK ,UK , VK , X J , Y J ,UJ , V J , X A , YA, UA , V A, DEL A, BET AK , S I GK , S IGSQ, 

3 X ( 50,50 ) , Y ( 50,50 ) ,U( 50,50) ,V ( 50,50 ) ,BET A ( 50 ) ,PSIR ( 50) , 

4 UI,V1 , DELI ,DELY,RPRES, RDENS,RECOV, ENT P,NDIM,KREF, I REF, 

5 XB (50) ,YB (50) »UB( 50) »VB( 50) ,P(50) ,RECV<2, 10) , FACTR , IREV , 

6 REG,R<25,4) ,S(25,4) , COWL ( 25 ) , BODY ( 25 ) , NR , NS , I COWL , I BODY 
DIMENSION NPT ( 2 ) 

K = KA 
I = I A 

IF (K ,E0. I) GO TO 30 
J F = K A- 1 A 
DO 10 J=1 , JF 
I = I A+ J 

IF ( X ( I , K ) .GT, 0.0) GO TO 12 
10 CONTINUE 
12 I REF= I 
XF = XR 
YF = YR 

DO 24 J = 1 ,2 
RM I N= 10.0 
DO 14 K= I REF ,K A 

IF ( X ( I R EF , K ) .EO. 0.0) GO TO 14 

RADI US = SQRT ( ( X ( I REF ,K ) -XF ) #*2+ ( Y ( I R EF ,K ) -YF ) **2 ) 

RMIN=AMIN1 (RADIUS, RMIN) 

14 CONTINUE 

DO 16 K=IREF,KA 

RADI US = SQRT ( ( X ( I REF ,K ) -XF ) **2 + < Y ( I R EF , K ) -YF ) ) 

IF (RADIUS .EO. RMIN) GO TO 18 
16 CONTINUE 
18 KREF=K 

IF (IREF .L T. KR EF > GO TO 20 
K=KREF+ 1 

IF ( X ( I REF ,K ) .EO. 0.0) GO TO 30 
GO TO 22 

20 CALL BOUND ( IREF, KREF, IREF, KREF-1 , XR , YR , SLOPE, SIGN) 

XREF=XK 

TEST = XREF-X ( I REF, KREF) 

KSIGN=TEST/ABS(TEST) 

K=KREF+KSIGN 

IF (IREF .GT. K ) GO TO 30 
IF ( X ( I R EF , K ) .EO. 0.0) GO TO 30 
22 CALL BOUND! IREF, KREF, IREF, K,XR,YR, SLOPE, SIGN) 

XF = XK 
YF = YK 

TEST=XF-X(KREF,IREF) 

KSIGN=TEST/ABS(TEST) 

K=KREF+KS IGN 

I F ( IREF .GT. K ) GO TO 30 
24 CONTINUE 

CALL FIND(XK) 

IF (REG .EO. -1.0) GO TO 26 
C0 = R ( I COWL , 1 ) 

C1=R( I COWL ,2 ) 

C2 = R ( I COWL , 3 ) 
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C3 = R ( ICOWL 1 4 ) 

GO TO 28 

26 CO = S( IBODY»l ) 

C1 = S( I BODY , 2 ) 

C2 = S ( I BODY , 3 ) 

C3= S ( I BODY , 4 ) 

28 YSURF=CO+Cl*XK+C2*XK**2+C3*XK**3 
TEST=( YK-YSURF )/ABS( YK-YSURF ) 

IF (TEST .EG. SIGN) GO Tn 56 
30 N P T ( 2 ) = 0 
XREF=XR 
DO 54 L= 1 , 2 
CALL FIND(XREF) 

IF (REG .EQ. -1*0) GO TO 32 
CO=R ( ICOWL, 1 ) 

C 1 = R ( ICOWL, 2 ) 

C2 = R ( ICOWL, 3) 

C3=R< ICOWL, 4) 

PO I NT = C OWL ( ICOWL) 

NPT ( L ) = ICOWL 
GO TO 34 
32 CO=S< I BODY , 1 ) 

Cl = S ( I BODY , 2 ) 

C2 = S ( I BODY , 3 ) 

C3 = S ( I BODY , 4 ) 

PO I NT=BODY ( I BODY ) 

NPT ( L ) = I BODY 

34 IF ( NPT ( 1 ) .EQ. NPT(2)) GO TO 56 
IF (ABS(SLOPE) .GT. ERROR) GO TO 36 
XK=XR 
GO TO 46 
36 AL PHA = C 3 
ETA=C2 
SIGA=C1 

DELTA=S IGA -SLOPE 

GAMMA= (CO-YR)+SLOPE*XR 

CALL R00T3 (ALPHA, ETA, DELTA, GAMMA) 

M = 0 

DO 38 1=1,3 
J=2*I~1 
K = 2# I 

P( J)=P(J)-POINT 

IF (P(J) .LE. 0.0) GO TO 38 

IF ( ABS ( P ( K ) ) .GT. ERROR) GO TO 38 

M=M+ 1 

PS I R ( M ) =P ( J ) 



38 CONTINUE 
XM I N= 1 0 , 0 
DO 40 1 = 1, M 

XMIN=AMIN1 ( XM I N, P.S I R ( I ) ) 

40 CONTINUE 
44 X K = X M I N+ POINT 

46 YK=C0+C1*XK+C2*XK**2+C3*XK**3 
SIGMA=C1+2,0*C2*XK+3,0*C3*XK**2 
RMIN=10.0 
DO 48 1=1,50 

IF ( XB ( I ) .EG. 0*0) GO TO 48 

TEST = SORT ( ( X B ( I ) -XK ) **2+ ( YB ( I )-YK)**2) 

RM I N = AM I N 1 (TEST t RM 1 N ) 

48 CONTINUE 

DO 50 1 = 1 t 5 0 

IF ( XB ( I ) • EG* 0,0) GO TO 50 
T EST = SORT ( ( XB ( I ) -XK ) **2+ < YB ( I >-YK)**2) 

IF (RMIN • EG . TEST) GO TO 52 
50 CONTINUE 
52 TEST = XK- XB ( I ) 

ISIGN=TEST/ABS( TEST ) 

J= 1+ I S I GN 
ANGLE=ATAN( S I GMA ) 

G I = SOR T ( UB ( I )**2+VB< I )** 2 ) 

0 J = SOR T ( UB ( J )**2+VB ( J)**2 ) 

DS=SORT( ( XB ( J )-XB< I ) ) **2+ ( YB ( J ) -YB < I ) )**2) 

DELS=SORT ( ( XK-XB ( I ) )**2+(YK-YB< I ))**2) 

DQ=QJ-0I 

DODS=DO/DS 

OK=QI+DODS*DELS 

UK=QK*COS< ANGLE) 

VK=OK*S IN( ANGLE ) 

KR EF = 50 
I R EF= I + 1 
XR EF = XK 
54 CONTINUE 
56 RETURN 
END 



$ I 8FTC LOCI LIST, REF, DECK, DEBUG 
C 

C SUBROUTINE FOR LOCATING A SPECIFIC CHARACTERISTIC 

C 

SUBROUTINE LOCUS < KA , S I GN ) 

COMMON AMO, GAM, THETC, DEL B , DELU1 , DELU, UC , VC, DELC,BS, DELS, ERROR, 

1 0,AM,A MU, THETA, A, B,C,D,E,F,G, PRES, DENS, AREA, DYNP, CP, PSI, 

2 XK,YK,UK,VK,XJ,YJ,UJ,VJ,XA,YA,UA,VA,DELA,BETAK,SIGK,SIGSQ, 

3 X( 50, 50),Y(50,50) ,U( 50,50) ,V( 50,50), BETA( 50) ,PSIR(50) , 

4 UI ,VI,DEL1,DELY,RPRES,RDENS,REC0V,ENTP,MDIM,KREF, I REF, 

5 XB ( 50 ) , YB ( 50 ) , UB ( 50 ) ,VB(50) , P ( 50 ) , R EC V ( 2 , 1 0 ) , F ACTR , I REV , 

6 REG,R(25,4) ,5(25,4) , COWL (25 ) , BODY ( 2 5 ) , NR , NS , I COWL, I BODY 
DO 8 K=K A , 50 

IF ( X(K,1 ) * EQ. 0.0) GO TO 10 
8 CONTINUE 
10 KREF=K 

DO 12 1=2,50 

IF ( X ( KREF-1 , I ) .GT. 0.0) GO TO 14 
12 CONTINUE 
14 I R EF= I 

KP=KREF-1 
REG=-REG 
XP = X ( KP , 1 ) +3 .0 
TANS=V(KP,1 )/U(KP,l) 

CALL CURVE ( X (KP, 1 ) , Y ( K P , 1 ) , TANS , XP , 0 .0 , 0 . 0 , 1 . 0 , NS ) 

NS=NS+1 

CALL CCRE(X(KP,IREF) ,Y(KP,IREF) ,U(KP,IREF),V(KP,IREF),-SIGN) 

CALL EOUK ( KREF , I REF , XK , YK , UK , VK ) 

CALL CGREtKREF, IREF, -1.0) 

REG=-REG 

BODY(NS )=X (KREF, IREF) 

CALL CFPR( X(KREF, IREF) ,-Y (KREF, IREF) ,U(KREF,'l REF) ,-V( KREF, IREF) ) 

PS I R ( I REF ) =PS I 

STFR= 1 . 0 

KP=KRE F-l 

DO 22 KK = 1 , K P 

K=KREF-KK 

IF ( IREF .GT. K) GO TO 24 

CALL CFPR(X(K,IREF) ,-Y (K, IREF) ,U( K , IREF ) ,-V ( K , IREF ) ) 

PSIAV=.5*( PSIR( IREF) +PS I ) 

PS I R ( IREF)=PSI 

IF ( NO I M .GT. 2) GO TO 18 

REF = Y ( K , I REF ) — Y ( K+l , IREF ) 

GO TO 20 

18 REF=Y(K,IREF)**2-Y(K+1, IREF)**2 
20 STFR=STFR+PSIAV*REF 
22 CONTINUE 
24 DO 28 K= 1 , 50 
DO 28 1=1,50 
IF ( I .GT. K ) GO TO 26 
IF (K .GT. KREF) GO TO 26 
IF ( I .LE. IREF) GO TO 28 
26 CALL EOUK (K ,1 ,0.0, 0.0, 0.0, 0.0) 

28 CONTINUE 

DO 30 1=1,50 

IF (I .LE. IREF) GO TO 30 
CALL EQUI ( I ,0.0, 0.0,0. 0,0.0) 

30 CONTINUE 
PS I =STFR 
RETURN 
tND 
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$ I BFTC SURFS L I S T , RE F , DECK , DEBUG 

CONDITIONS IN THE I S5NTROP I C COMPRESSION REGION 

SUBROUTINE SURF (KS, I S , NN , AMT , THROAT , ANGLE , STFR , S I GN ) 

COMMON AMO , GAM , THETC, DEL B, DELU 1,DELU,UC, VC, DELC,BS, DELS, ERROR, 

1 0,AM,AMU,THETA,A,B,C,D,E,F,G,PRES.,DENS,AREA, DYNP,CP,PSI, 

2 XK , YK ,UK,VK,XJ,YJ,UJ,VJ ,XA,YA,UA,VA,DELA,BETAK, SIGK,SIGSQ, 

3 X( 50,50) ,Y< 50,50) ,U< 50,50) ,V< 50,50) ,BETA( 50) ,PSIR( 50) , 

4 U1 ,V1 ,OELI ,OELY,RPRES,RDENS,RECOV,ENTP,MDIM,KREF, I REF, 

5 XB( 50) , YB( 50) ,UB( 50) ,VB( 50) , P ( 50 ) , R EC V ( 2 , 1 0 ) , F ACTR , I RE V , 

6 REG,R (25,4) , S ( 25 ,4 ) , COWL ( 25 ) , BODY ( 25 ) , NR , NS , ICOWL , I BODY 
100 FORMAT ( / / 8X , 1 HX , 1 3X , 1 HY , 1 3X , 1H0, 1 1 X , 5HTHET A, 1 1 X , 1 HU, 1 3X , 1HV , 9X , 

18HMACH NO. ,8X,4HP/H0, 10X,4HP/P0// ) 

102 FORMAT (1P9E14.5//) 

104 FORMAT ( / /43X , 33HC OND I T IONS ON THE CONTOUR SURFACE//) 

106 FORMAT ( //45X,28HC0NDI T IONS IN THE FLOW FIELD//) 

108 FORMAT ( / / 44X , 3HAMT , 1 1 X , 3HTHR , 1 1 X , 3HANG/ / ) 

110 FORMAT ( 37X,1P3E14.5// ) 

112 FORMAT ( //28X,46HISENTR0PIC COMPRESSION REGION INITIATES AT X =1PE 
112.5, 2X, 3HY =1 PE 12. 5////) 

114 FORMAT ( 1 H 1 ) 

116 FORMAT (1HJ) 

118 FORMAT <//33X,54HCALCULATI0N OF FLOW FIELD PERFORMED IN SUBROUTINE 
1 SURF//) 

120 FORMAT ( / /29X , 50HPRESSURE RECOVERY (H/HO) ALONG THE UPPER CONTOUR 
1 = 1 PE 12. 5//) 

122 FORMAT (/ /29X , 50HPRESSURE , RECOVERY (H/HO) ALONG THE LOWER CONTOUR 
1=1PE12.5// ) 

K = KS 
I = IS 
MP = K 

MPF=MP+NN 
XNN=NN 
PO I NT= 1 .0 

CONVR= ( 1 .0+ ( GAM-1.0 )/2.0*AM0**2 )**( -GAM/ ( GAM-1.0) ) 

I REV= I REV+1 

RECV ( 1 , IREV )=RECV ( 1 , IREV-1 ) 

RECVI2, IREV)=RECV(2, 1 REV-1 ) 

CALL CFPR(X(K.I),Y(K,I),U(K,I),V(K,I)) 

AMR = AM 

THETAR=THETA 
SIGI=TAN(THETAR) 

IP=I+1 
P( I ) = 1 . 0 

DELAM=( AMR-AMT ) / XNN 
DELX=THROAT/XNN 
TAU=TAN( ANGLE ) -TAN ( THE TAR ) 

XO=X(KS, IS) 

AO=Y(KS,IS) 

A1 = V(KS» IS)/U(KS, IS) 

IF (THROAT .EQ. 0.0) GO TO 10 
A2=TAU/(2.0*THR0AT) 

10 AMF=AMR 

NNN=NN+I P-1 
DO 20 I = I P , NNN 

IF (THROAT .GT. 0.0) GO TO 12 
P ( I ) = 1 .0 
GO TO 14 

12 P ( I ) =1 .0+P ( 1-1 ) 
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14 AMF=AMF-DELAM 
KP = P ( I ) 

KT=P( I ) — P ( 1-1 ) 

K=MP— 1+KP 
KK=K-KT 
KL = K 

XK = X ( KK , 1-1 ) +OEL X 

IF (ABS(ANGLE) .GT. 0.0) GO TO 16 
CALL PMER(AMR,-THETAR,AMF) 

S I GK=— VK/UK 
SIGC=0.5*(SIGI+SIGK) 

S I G I =S I GK 

YK=Y (KK, 1-1 )+SIGC*DELX 

EXP=3.0 

GO TO 18 

16 YK=A0+A1*( XK-X0)+A2*( XK-X0)*#2 
SIGK=A1+2.0*A2*( XK-XO) 

PHI =AT AN ( S IGK ) 

OK=SORT ( AMF**2/ ( 1 . 0+S I GSO* ( AMF**2-1 . 0 ) ) ) 

UK=QK*COS(PHI ) 

VK=-QK*SIN< PHI ) 

EXP = 2 .0 

18 CALL EOUMK, I , XK , YK , UK ,-VK ) 

KREF=K 
I REF= I 

CALL CFPRIXIK, I ) ,-Y(K,I ) ,U(K, I ) ,-VIK, I ) ) 

PS I R ( I ) = PS I 
20 CONTINUE 

DO 36 I = I P , NNN 
STF=1 .0 

DO 34 KK=1,MPF 
KP = P ( I )-l .0 
K = MP+ 1 +K P— KK 

22 CALL CCRCI X (K— 1 , I— 1 ) ,Y(K-1,I-1 ) ,U(K-1 , 1-1 ) , V(K— 1, 1-1 ) ,X (K r I ) t 
1 Y(K, I ) ,U(K, I ) ,V< K, I ) ,SIGN) 

CALL E0UK(K-1,I ,XK,YK,UK,VK) 

24 CALL CFPR ( X(K-1 , I ) ,-Y(K-l, I ) ,U(K-1, I ) ,-V(K-l, I ) ) 

PS I A V= . 5* ( PS I R ( I ) + PSI ) 

IF ( NO I M .GT. 2) GO TO 26 
REF = Y(K — 1,1 ) — Y ( K , I ) 

GO TO 28 

26 REF=Y(K-1 ,1 )**2-Y(K,I )*%2 
28 STF=STF+PSIAV*REF 
PSIRI I )=PSI 

IF (STFR-STF) 30,32,34 

30 CALL CCSRI STFR , STF , X ( K-l , I ) , -Y ( K-l , I ) , U ( K-l , I ) ,-V ( K-l , I ) , X ( K , I ) , 
1 -Y(K, I ) ,U(K, I ) ,-V(K, I ) ) 

CALL EOUII I ,XK,-YK,UK,-VK) 

GO TO 36 

32 CALL EOUI ( I ,X(K-1 ,1 ) ,Y(K— 1, I ) ,U(K-1 ,1 ) ,V (K— 1 ,1 ) ) 

GO TO 36 
34 CONTINUE 
36 CONTINUE 

DO 38 J = 1 ,KS 
K= ( KS+ I )-J 

IF ( X ( K » I S ) .EO. 0.0) GO TO 40 
38 CONTINUE 
40 KR = K+ 1 
I R= I S 
R EG = -1 . 0 
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CALL CURVE ( X(KS, IS) ,Y(KS, IS) , V ( KS, I S ) /U( KS, I S ) , X ( KREF, I REF) , 

1 Y ( KRE F , I R E F ) , V ( KREF , 1 REF ) /U ( KREF , I REF ) , EXP , NS ) 

NS=NS+ 1 
R EG= 1 • 0 

CALL CURVE (X6( IS ) ,YB( IS) ,VB< IS)/UB( IS),XB< IREF) , YB(IREF) , 

1 VB( IREF)/UB( I REF) ,3.0, NR) 

NR=NR+ 1 
WRITE (6,114) 

WRITE (6,118) 

WRITE (6,108) 

WRITE (6,110) AMT, THROAT, ANGLE 
WRITE (6,120) R EC V ( 2 , I R E V ) 

WRITE (6,122) RECV(1,IREV) 

WRITE (6,104) 

WRITE (6,100) 

DO 42 I =1 R , I REF 

IF (XB( I ) .EO. 0.0) GO TO 42 

CALL CFPR(XB( I ) ,YB( I ) ,UB( I ) ,VB( I ) ) 

PRHO=FACTR*PRES 

PRPO=PRHO/CGNVR 

WRITE (6,102) X B ( I ),YB( I ),0, THETA, UB( I),VB( I ) , AM , PRHO , PR PO 
42 CONTINUE 

WRITE (6,114) 

WRITE (6,106) 

WRITE (6,100) 

DO 46 K=KR,KR6F 
WRITE (6,116) 

DO 46 I=IR, IREF 

IF ( X ( K , I ) .EO. 0.0) GO TO 46 

T E ST = AB S ( X ( K , I ) —X ( K S , I S ) ) 

IF (TEST .GT. 0.0) GO TO 44 
IF (POINT .EO. 0.0) GO TO 44 
WRITE (6,112) X(K, I ) ,Y(K, I ) 

POINT = 0 .0 

44 CALL CFPR( X(K, I ) , Y (K, I ) ,U(K, I ) ,V ( K, I ) ) 

PRHO=FACTR*PRES 

PRPO=PRHO/CONVR 

WRITE (6,102) X (K, I ) ,Y(K,I ) , 0 , THETA , U ( K , I ) ,V(K, I ) , AM , PRHO , PR PO 

46 CONTINUE 

IF ( XB ( IREF) *LT*X(KS,IS) ) GO TO 47 
CALL ENDI (KS, IS, SIGN) 

47 DO 48 J = 1 ,KREF 
K= ( KREF+ 1 )~J 

IF ( X ( K , I R E F ) .EO. 0.0) GO TO 50 

48 CONTINUE 
50 KP=K+ 1 

DO 52 J= 1 , 50 
K= ( K P — 1 )+J 

IF ( X(K, IREF) .EQ. 0.0) GO TO 54 

CALL EOUK( J ,1 ,X(K, IREF ),Y(K, IREF), U ( K , I REF ) , V ( K , IREF)) 

52 CONTINUE 
54 KREF=J^1 

CALL E OU I ( 1 , X B ( IREF) , YB ( I R EF ) , UB ( I RE F ) , VB ( I R E F ) ) 

I REF=1 

DO 58 K= 1 , 50 
DO 58 1=1,50 

IF (K .GT. KREF) GO TO 56 
IF (I .EQ. 1 ) GO TO 58 
56 CALL E0UK(K,I,0.0,0.0,0.0,0.0) 

58 CONTINUE 

OU 60 1=2,50 

CALL EOUI ( I ,0.0, 0.0, 0.0, 0.0) 

60 CONTINUE 
RETURN 
END 
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$ I BFTC ENDIR L I ST , REF , DECK , DEBUG 

C 

C 

SUBROUTINE END I ( K A , I. A * S I GN ) 

COMMON AMO,GAM,THETC,DELB,DELUl ,DELU,UC , VC, DELC,BS, DELS , ERROR t 

1 0, AM, AMU, THETA, A, B,C, D, E, F ,G,PRES, DENS , AREA, DYNP,CP, PS I , 

2 XK , YK ,IJK , VK , X J , V J , U J , V J ,X A , Y A , UA , VA , DEL A , BETAK , S I GK , S I GSQ , 

3 X(50,50),Y(50,50) , U ( 50 , 50 ) ,V(50,50),BETA(50),PSIR(50) , 

4 U1 ,V1 , DEL 1, DEL Y, RPRES, R DENS, R ECOV, ENT P,NDIM,KREF, I REF, 

5 XB( 50) ,YB( 50) ,UB( 50) ,VB(50) ,P ( 50) ,RECV(2, 10) ,FACTR, I REV, 

6 REG,R(25,4) ,S( 25,4) , COWL ( 25) ,BODY(25) , NR T NS , ICOWL , I BODY 
100 FORMAT (//38X,46HC OND I T I ONS IN THE VICINITY OF THE NORMAL SHOCK//) 
102 FORMAT { / / 29X , 1 HX , 1 2X , 1 H Y , 1 OX , 8HMACH NO . , 8 X , 5HTHET A , 9 X , 5HP2 / P 1 , 9X , 

14HH/H0// ) 

104 FORMAT ( 2 1 X , 1 P6E 14 . 5/ / ) 

200 FORMAT ( 1 HI ) 

P I =3 • 1 41 5927 

CONVP=< 1.0+ (GAM-1. 0) /2.0)** (-GAM/ (GAM- 1.0) ) 

RAVE=FACTR/CONVP 
KR=K A- 1 

CALL CFPR(X(KA,IA) ,Y (KA, I A) ,U(KA, IA) ,V(KA, I A) ) 

SLOPF=TAN( PI/2* 0+ THETA) 

CALL OSWRUM, PI/2.0) 

RECOV=RAVE*RECOV 
WRITE (6,200) 

WRITE (6,100) 

WRITE (6,102) 

WRITE (6,104) X(KA, IA) , Y ( K A , I A ) , AM , THE T A , R PR ES , R ECOV 
XR E F = X ( K A , I A ) 

YREF=Y ( K A , I A ) 

DO 18 KK= 1 , KR 

K=K A-KK 

XR = XRE F 

YR = YRF F 

DO 16 J = 1 , 2 

RMIN=10.0 

DO 10 I = 1 , I REF 

IF ( X ( K , I ) .EO. 0.0) GO TO 10 

RADI US = SOR T ( ( X ( K , I ) -XR ) **2+ ( Y ( K , I )-YR)**2) 

RMIN=AMIN1(RM IN, RADIUS) 

10 CONTINUE 

DO 12 1 = 1, I REF 

IF ( X ( K , I ) .EO. 0.0) GO TO 12 

RADI US=SORT ( (X (K, I ) -XR ) **2+ ( Y ( K , I)-YR)**2 ) 

IF (RMIN .EO. RADIUS) GO TO 14 
12 CONTINUE 

14 IF ( X ( K , I + 1 ) .EO. 0.0) GO TO 20 

CALL BOUND (K, I , K , I + 1 , XRE F , YR EF , SLOPE , S I GN ) 

TEST=XK-X(K, I ) 

ISIGN=TEST/ABS(TEST) 

I P= I + I S IGN 

IF ( X(K, IP) .EO. 0.0) GO TO 20 

CALL BOUND (K, I ,K , I P , XREF , YREF , SLOPE , S I GN ) 

X R = XK 
YR = YK 

16 CONTINUE 
XR EF = XR 
YR E F = YR 

CALL CFPR( XK , YK ,UK , VK ) 

SLOPE=TAN( PI/2.0+THETA) 
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CALL OSWRt AM, PI/2.0) 

RECOV=RAVE*RECOV 

WRITE (6,104) XK,YK, AM, THETA, RPRES,RECOV 
18 CONTINUE 
20 RMIN=10.0 

DO 22 I = 1 , I REF 

R AD I US = SQR T ( ( XB ( I > -XREF >**2+( YB < I)-YREF)**2) 

RM IN* AM INI ( RM IN, RAD I US ) 

22 CONTINUE 

DO 24 1 = 1, I R E F 

RADI US = SOR T ( ( XB ( I ) -XREF ) **2+ ( Y B ( I )-YREF)**2) 

IF (RMIN .EO. RADIUS) GO TO 26 
24 CONTINUE 
26 TEST=XREF-XB(I ) 

ISIGN=TEST/ABS(TEST ) 

J = I + I S I GN 
DU=UB ( J ) —UB ( I ) 

DV = VB ( J ) — VB ( I ) 

DX = XB( J ) — XB ( I ) 

DY=YB( J)-YB( I ) 

DS=SORT ( DX**2+DY**2 ) 

DUOS=DU/OS 
DYDX=DY / DX 

XK= ( YB ( I )-YREF+SLOPE*XREF-DYDX*XB( I ))/( SLOPE-O YDX ) 
YK= YB ( I )+DYDX#(XK-XB{ I ) ) 

DELS=SORT( ( XK-XB ( I ) )**2+(YK-YB< I ) )** 2 ) 

UK =UB ( I )+DUDS*DELS 
VK=UK*DYDX 

CALL CFPR( XK,YK,UK,VK) 

SLDPE=TAN( PI/2.0+THETA) 

CALL OSWR( AM, PI/2.0) 

RECOV=RAVE*RECOV 

WRITE (6,104) XK,YK, AM, THETA, RPRES,RECOV 

RETURN 

END 



$ I BF TC ENDSR L I ST, REF, DECK, DEBUG 

C 

C 

SUBROUTINE ENDS ( K A , I A , S I GN ) 

COMMON AMO, GAM, THETC, DELB, DELU1 , DELU, UC , VC, DELC,BS, DELS, ERROR, 

1 O, AM, AMU, THETA, A, B,C,D,E,F,G, PRES, DENS, AREA, DYNP, CP, PSI, 

2 XK,YK,UK,VK,XJ, YJ,UJ, V J , X A , Y A , UA , V A , DEL A , BE T AK , S I GK , S I GSO , 

3 X( 50,50) ,Y<50, 50) ,U( 50,30) , V ( 50,50 ) , BETA (50) ,PSIR( 50) , 

4 U1 , VI, DELI ,DELY,RPRES,RDENS,RECOV,ENTP,NDIM,KREF, I REF, 

5 XB(50) ,YB(50),UB(50) , VB ( 50 ) , P ( 50 ) ,RECV(2,10) , F ACTR , I REV , 

6 REG,R(25,4) , S ( 2 5 , 4 ) , COWL (25) , BODY ( 2 5 ) , NR , NS , I COWL , I BODY 
100 FORMAT ( / / 38X , 46HC0ND I T I ONS IN THE VICINITY OF THE NORMAL SHOCK//) 
102 FORMAT ( / /2 9X , 1 HX , 1 2 X , 1 HY , 1 OX , 8HM ACH NO . , B X , 5HTHET A , 9X , 5HP2 / P 1 , 9X , 

14HH/H0// > 

104 FORMAT ( 2 1 X , 1 P6E 1 4 . 5 / / ) 

200 FORMAT ( 1 HI ) 

TEST=0. 01745329 
P I =3 • 1415927 

C ONVP = ( 1.0+ (GAM-1'. 0>/2.0)** (-GAM/ (GAM- 1.0) ) 

RAVE = FACT R /C ON VP 
KREF=K A 
I REF = I A 

XK = X ( I REF , KREF ) 

YK = Y ( IREF,KREF ) 

UK = U ( I R EF , KREF ) 

VK = V( IREF,KREF) 

X P = XK 
YP = YK 

CALL CFPR(XK,YK, UK , VK ) 

IF (ABS(THETA) .GT. TEST) GO TO 6 
SL0PE=0 • 0 
GO TO 8 

6 SLOPE=TAN( PI/2.0+THETA) 

8 CALL OS WR (AM, PI/2.0) 

RECOV=RAVE*RECOV 
WRITE (6,200) 

WRITE (6,100) 

WRITE (6,102) 

WRITE (6,104) XK,YK, AM, THETA, RPRES,RECOV 
I F=K A —I 
DO 14 1 = 1, IF 

IF (KREF .EQ. 50) GO TO 16 

CALL SEEK (KREF, I REF, XP,YP, SLOPE, SIGN) 

CALL CFPR(XK,YK,UK,VK) 

IF ( ABS ( THE TA ) .GT. TEST) GO TO 10 
S L0PE = 0 .0 
GO TO 12 

10 SLOPE=TAN( PI/2.0+THETA) 

12 CALL OSWR(AM, PI/2.0) 

R ECOV = R A V E^R ECOV 

WRITE (6,104) XK,YK, AM, THETA, RPRES,RECOV 
X P= XK 
Y P = YK 

14 CONTINUE 
16 RETURN 
END 
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Initial 
datum line 


C_ characteristic 



Figure 1. - Domain in the x,y plane in which solution can 
be established. 



Figure 2. - Domain in x, y plane in which solution can be established 
for free-boundary problem. 





x-coordinate 

Figure 3. - Location of basic field point. 



x(K-U-l), y(K-U-l) 









Figure 7. - Location of shock-boundary points. 


Cowl surface, y = y c (x) 



x-coordinate 

Figure 8. - Fixed-boundary problem. 




x-coordfnate 

Figure 9. - Free-boundary problem. 
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Initial characteristic 





o Centerbody surface 
□ Cowl surface 



(a) Mach number distribution. 
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Dimensionless x-coordinate 
(b) Characteristic solution. 

Figure 11. - Bicone mixed-compression inlet (10.0°-18.5°) designed for free-stream Mach number of 2.500. Cowl angle, 5.0°; throat Mach 
number, 1.300. 







Dimensionless y-coordinate Mach number 






Dimensionless y-coordinate Static-pressure ratio. p/P Q 





(c) Schlieren photograph showing isentropic spike cone. 


Figure 15. - Isentropic spike configuration designed for free-stream Mach number of 2. 49 and focal point Mach number of 1 46 Initial cone 
angle, 16. 0 . 
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Dimension less y-coordinate Static-pressure ratio, p/ P q 







to) Schlieren photograph of bicone configuration. 

Figure 17. - Bicone forebody (10°-18.5°) configuration for free-stream Mach number of 2.49. 
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